TABLE OF CONTENTS
ABINIT/chi0_bksmask [ Functions ]
NAME
chi0_bksmask
FUNCTION
Compute tables for the distribution and the storage of the wavefunctions in the SCREENING code.
INPUTS
Dtset<type(dataset_type)>=all input variables for this dataset Ep<em1params_t>=Parameters for the screening calculation. Kmesh <kmesh_t>=Structure describing the k-point sampling. nbvw = Max. number of fully/partially occupied states over spin nbcw = Max. number of unoccupied states considering the spin nprocs=Total number of MPI processors my_rank=Rank of this this processor.
OUTPUT
bks_mask(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state. keep_ur(Ep%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space. ierr=Exit status.
PARENTS
screening
CHILDREN
xmpi_split_work2_i4b
SOURCE
718 subroutine chi0_bksmask(Dtset,Ep,Kmesh,nbvw,nbcw,my_rank,nprocs,bks_mask,keep_ur,ierr) 719 720 use defs_basis 721 use defs_datatypes 722 use defs_abitypes 723 use m_gwdefs 724 use m_errors 725 use m_xmpi 726 727 use m_bz_mesh, only : kmesh_t 728 729 !This section has been created automatically by the script Abilint (TD). 730 !Do not modify the following lines by hand. 731 #undef ABI_FUNC 732 #define ABI_FUNC 'chi0_bksmask' 733 !End of the abilint section 734 735 implicit none 736 737 !Arguments ------------------------------------ 738 !scalars 739 integer,intent(in) :: my_rank,nprocs,nbvw,nbcw 740 integer,intent(out) :: ierr 741 type(Dataset_type),intent(in) :: Dtset 742 type(em1params_t),intent(in) :: Ep 743 type(kmesh_t),intent(in) :: Kmesh 744 !arrays 745 logical,intent(out) :: bks_mask(Ep%nbnds,Kmesh%nibz,Dtset%nsppol) 746 logical,intent(out) :: keep_ur(Ep%nbnds,Kmesh%nibz,Dtset%nsppol) 747 748 !Local variables------------------------------- 749 !scalars 750 integer :: my_nspins,my_maxb,my_minb,isp,spin,distr_err,nsppol,band,rank_spin,ib 751 character(len=500) :: msg 752 logical :: store_ur 753 !arrays 754 integer :: my_spins(Dtset%nsppol),nprocs_spin(Dtset%nsppol) 755 integer,allocatable :: istart(:),istop(:) 756 757 ! ************************************************************************* 758 759 ierr=0; nsppol=Dtset%nsppol 760 761 my_nspins=Dtset%nsppol; my_spins= [(isp,isp=1,nsppol)] 762 763 ! List of spins for each node, number of processors per each spin 764 ! and the MPI rank in the "spin" communicator. 765 nprocs_spin = nprocs; rank_spin = my_rank 766 767 if (nsppol==2.and.nprocs>1) then 768 ! Distribute spins (optimal distribution if nprocs is even) 769 nprocs_spin(1) = nprocs/2 770 nprocs_spin(2) = nprocs - nprocs/2 771 my_nspins=1; my_spins(1)=1 772 if (my_rank+1>nprocs/2) then 773 my_spins(1)=2 774 rank_spin = my_rank - nprocs/2 775 end if 776 end if 777 778 store_ur = (MODULO(Dtset%gwmem,10)==1) 779 bks_mask=.FALSE.; keep_ur=.FALSE. 780 781 select case (Dtset%gwpara) 782 case (1) 783 ! Parallelization over transitions **without** memory distributions (Except for the spin). 784 my_minb=1; my_maxb=Ep%nbnds 785 do isp=1,my_nspins 786 spin = my_spins(isp) 787 bks_mask(my_minb:my_maxb,:,spin) = .TRUE. 788 if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE. 789 end do 790 791 case (2) 792 ! Distribute bands and spin. 793 do isp=1,my_nspins 794 spin = my_spins(isp) 795 796 if (nprocs_spin(spin) <= nbcw) then 797 ! Distribute nbcw empty bands among nprocs_spin (block of bands without replicas). 798 ! Bands are distributed in contiguous blocks because 799 ! this distribution is well suited for the Hilber transform 800 ! since each node will allocate only a smaller frequency interval 801 ! for the spectral function whose size scales with the number of MPI nodes. 802 ! Note it is now meaningless to distinguish gwcomp=0 or 1 since the workload is well balanced later on 803 ABI_MALLOC(istart,(nprocs_spin(spin))) 804 ABI_MALLOC(istop,(nprocs_spin(spin))) 805 806 call xmpi_split_work2_i4b(nbcw,nprocs_spin(spin),istart,istop,msg,distr_err) 807 808 if (distr_err==2) then 809 ! Too many processors. 810 !MSG_WARNING(msg) 811 ierr=1 812 end if 813 814 my_minb = nbvw + istart(rank_spin+1) 815 my_maxb = nbvw + istop (rank_spin+1) 816 817 ABI_FREE(istart) 818 ABI_FREE(istop) 819 820 if (my_maxb-my_minb+1<=0) then 821 write(msg,'(3a,2(i0,a),2a)')& 822 & 'One or more processors has zero number of bands ',ch10,& 823 & 'my_minb= ',my_minb,' my_maxb= ',my_maxb,ch10,& 824 & 'This is a waste, decrease the number of processors.' 825 MSG_ERROR(msg) 826 end if 827 828 bks_mask(my_minb:my_maxb,:,spin)=.TRUE. 829 if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE. 830 831 else 832 ! New version (alternate bands with replicas if nprocs > nbcw) 833 ! FIXME: Fix segmentation fault with Hilbert transform. 834 do ib=1,nbcw 835 if (xmpi_distrib_with_replicas(ib,nbcw,rank_spin,nprocs_spin(spin))) then 836 band = ib + nbvw 837 bks_mask(band,:,spin)=.TRUE. 838 if (store_ur) keep_ur(band,:,spin)=.TRUE. 839 end if 840 end do 841 end if 842 843 ! This is needed to have all the occupied states on each node. 844 bks_mask(1:nbvw,:,spin) = .TRUE. 845 if (store_ur) keep_ur(1:nbvw,:,spin)=.TRUE. 846 end do ! isp 847 848 case default 849 ierr = 1 850 MSG_WARNING("Wrong value for gwpara") 851 end select 852 853 end subroutine chi0_bksmask
ABINIT/setup_screening [ Functions ]
NAME
setup_screening
FUNCTION
Initialize the Ep% data type containing the parameters used during the screening calculation. as well as basic objects describing the BZ sampling .... TODO list to be completed
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
wfk_fname=Name of the input WFK file. acell(3)=length scales of primitive translations (Bohr). rprim(3,3)=dimensionless real space primitive translations. ngfftf(18)=Contain all needed information about the 3D FFT for densities and potentials. dtfil <type(datafiles_type)>=variables related to files
OUTPUT
ngfft_gw(18)=Contain all needed information about the 3D FFT for the oscillator strengths. See ~abinit/doc/variables/vargs.htm#ngfft Ltg_q(:)<littlegroup_t>,= Ep<em1params_t>=Parameters for the screening calculation. Most part of it is Initialized and checked. Hdr_wfk type(Hdr_type)=Header of the KSS file. Cryst<crystal_t>=Definition of the unit cell and its symmetries. Kmesh<kmesh_t>=Structure defining the k-point sampling (wavefunctions). Qmesh<kmesh_t>=Structure defining the q-point sampling (screening) Gsph_wfn<gsphere_t>=Structure defining the G-sphere for the wavefunctions (not k-dependent). Gsph_epsG0<gsphere_t>=The G-sphere for the screening, enlarged to take into account for umklapps. Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file Vcp <type vcoul_t> datatype gathering information on the coulombian cutoff technique comm=MPI communicator.
SIDE EFFECTS
Dtset<Dataset_type>=All input variables for this dataset. %ecutwfn, %npwwfn, %ecuteps, %npweps might be redefined in setshells in order to close the shell.
PARENTS
screening
CHILDREN
xmpi_split_work2_i4b
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine setup_screening(codvsn,acell,rprim,ngfftf,wfk_fname,dtfil,Dtset,Psps,Pawtab,& 60 & ngfft_gw,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,KS_BSt,Ltg_q,Gsph_epsG0,Gsph_wfn,Vcp,Ep,comm) 61 62 use defs_basis 63 use defs_datatypes 64 use defs_abitypes 65 use defs_wvltypes 66 use m_errors 67 use m_profiling_abi 68 use m_xmpi 69 use m_nctk 70 #ifdef HAVE_NETCDF 71 use netcdf 72 #endif 73 use m_hdr 74 75 use m_gwdefs, only : GW_TOLQ0, GW_TOLQ, GW_Q0_DEFAULT, czero_gw, em1params_t 76 use m_fstrings, only : strcat, sjoin, ltoa, itoa 77 use m_io_tools, only : file_exists 78 use m_geometry, only : normv, mkrdim, metric 79 use m_crystal, only : crystal_print, crystal_t 80 use m_crystal_io, only : crystal_from_hdr 81 use m_bz_mesh, only : kmesh_t, kmesh_init, get_ng0sh, kmesh_print, find_qmesh, get_BZ_item,& 82 & littlegroup_t, littlegroup_init, make_mesh, kmesh_free 83 use m_ebands, only : ebands_init 84 use m_vcoul, only : vcoul_t, vcoul_init 85 use m_fft_mesh, only : setmesh 86 use m_gsphere, only : gsphere_t, gsph_init, merge_and_sort_kg, gsph_free, setshells 87 use m_pawtab, only : pawtab_type 88 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free 89 use m_io_kss, only : make_gvec_kss 90 use m_wfk, only : wfk_read_eigenvalues 91 92 !This section has been created automatically by the script Abilint (TD). 93 !Do not modify the following lines by hand. 94 #undef ABI_FUNC 95 #define ABI_FUNC 'setup_screening' 96 use interfaces_14_hidewrite 97 use interfaces_56_io_mpi 98 !End of the abilint section 99 100 implicit none 101 102 !Arguments ------------------------------------ 103 !scalars 104 integer,intent(in) :: comm 105 character(len=6),intent(in) :: codvsn 106 character(len=fnlen),intent(in) :: wfk_fname 107 type(Dataset_type),intent(inout) :: Dtset !INOUT is due to setshells 108 type(datafiles_type),intent(in) :: dtfil 109 type(Pseudopotential_type),intent(in) :: Psps 110 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw) 111 type(em1params_t),intent(out) :: Ep 112 type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out 113 type(ebands_t),intent(out) :: KS_BSt 114 type(kmesh_t),intent(out) :: Kmesh,Qmesh 115 type(crystal_t),intent(out) :: Cryst 116 type(gsphere_t),intent(out) :: Gsph_epsG0,Gsph_wfn 117 type(vcoul_t),intent(out) :: Vcp 118 !arrays 119 integer,intent(in) :: ngfftf(18) 120 integer,intent(out) :: ngfft_gw(18) 121 real(dp),intent(in) :: acell(3),rprim(3,3) 122 type(littlegroup_t),pointer :: Ltg_q(:) 123 124 !Local variables------------------------------- 125 !scalars 126 integer,parameter :: NOMEGAGAUSS=30,NOMEGAREAL=201,pertcase0=0,master=0 127 integer :: bantot,ib,ibtot,ikibz,iq,iqp,isppol,ig,ng,ierr 128 integer :: jj,mod10,mband,ng_kss,iqbz,isym,iq_ibz,itim 129 integer :: timrev,use_umklp,ncerr 130 integer :: npwepG0,nshepspG0,method,enforce_sym,nfftgw_tot !,spin,band,ik_ibz, 131 integer :: istart,iend,test_npwkss,my_rank,nprocs !ii 132 real(dp),parameter :: OMEGAERMAX=100.0/Ha_eV 133 real(dp) :: ecutepspG0,ucvol,domegareal 134 logical :: remove_inv,ltest,found,is_static,has_q0 135 character(len=500) :: msg 136 type(wvl_internal_type) :: wvl 137 !arrays 138 integer :: ng0sh_opt(3) 139 integer,allocatable :: npwarr(:) 140 integer,pointer :: gvec_kss(:,:) 141 integer,pointer :: test_gvec_kss(:,:) 142 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),qtmp(3),sq(3),qbz(3) 143 real(dp),pointer :: energies_p(:,:,:) 144 real(dp),allocatable :: doccde(:),eigen(:),occfact(:) 145 type(Pawrhoij_type),allocatable :: Pawrhoij(:) 146 147 ! ************************************************************************* 148 149 DBG_ENTER('COLL') 150 151 ! Check for calculations that are not implemented 152 ltest = ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol) == Dtset%nband(1)) 153 ABI_CHECK(ltest, 'dtset%nband(:) must be constant in the GW code.') 154 155 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 156 157 call mkrdim(acell,rprim,rprimd) 158 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 159 160 ! Set up basic parameters of the calculation 161 Ep%gwcalctyp =Dtset%gwcalctyp 162 Ep%plasmon_pole_model =.TRUE. 163 Ep%analytic_continuation=.FALSE. 164 Ep%contour_deformation =.FALSE. 165 166 mod10=MOD(Ep%gwcalctyp,10) 167 if (mod10/=0.and.mod10/=8) Ep%plasmon_pole_model =.FALSE. 168 if (mod10==1) Ep%analytic_continuation=.TRUE. 169 if (mod10==2.or.mod10==9) Ep%contour_deformation =.TRUE. 170 is_static=(mod10==5.or.mod10==6.or.mod10==7) 171 172 Ep%nbnds =Dtset%nband(1) 173 Ep%symchi =Dtset%symchi 174 Ep%inclvkb=Dtset%inclvkb; if (Dtset%usepaw/=0) Ep%inclvkb=0 175 Ep%zcut =Dtset%zcut 176 177 write(msg,'(2a,i4,2a,f10.6,a)')ch10,& 178 & ' GW calculation type = ',Ep%gwcalctyp,ch10,& 179 & ' zcut to avoid poles in chi0 [eV] = ',Ep%zcut*Ha_eV,ch10 180 call wrtout(std_out,msg,'COLL') 181 182 Ep%awtr =Dtset%awtr 183 Ep%npwe =Dtset%npweps 184 Ep%npwwfn=Dtset%npwwfn 185 Ep%npwvec=MAX(Ep%npwe,Ep%npwwfn) 186 187 timrev = 2 ! This information is not reported in the header 188 ! 1 --> do not use time-reversal symmetry 189 ! 2 --> take advantage of time-reversal symmetry 190 if (any(dtset%kptopt == [3, 4])) timrev = 1 191 192 if (timrev==1.and.Dtset%awtr/=0) then 193 MSG_ERROR("awtr/=0 cannot be used when time-reversal symmetry doesn't hold") 194 end if 195 196 ! Read parameters from WFK and verifify them. 197 call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm) 198 mband = MAXVAL(Hdr_wfk%nband) 199 call hdr_vs_dtset(Hdr_wfk,Dtset) 200 remove_inv=.FALSE. 201 202 test_npwkss = 0 203 call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,& 204 & gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr) 205 ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss") 206 207 ABI_MALLOC(gvec_kss,(3,test_npwkss)) 208 gvec_kss = test_gvec_kss 209 ng_kss = test_npwkss 210 211 if (Ep%npwvec>ng_kss) then 212 Ep%npwvec=ng_kss 213 if (Ep%npwwfn> ng_kss) Ep%npwwfn=ng_kss 214 if (Ep%npwe > ng_kss) Ep%npwe =ng_kss 215 write(msg,'(3a,3(a,i6,a))')ch10,& 216 & ' Number of G-vectors found less then required. Calculation will proceed with ',ch10,& 217 & ' npwvec = ',Ep%npwvec,ch10,& 218 & ' npweps = ',Ep%npwe ,ch10,& 219 & ' npwwfn = ',Ep%npwwfn,ch10 220 MSG_WARNING(msg) 221 end if 222 223 ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2)) 224 ierr = 0 225 do ig=1,ng 226 if (ANY(gvec_kss(:,ig)/=test_gvec_kss(:,ig))) then 227 ierr=ierr+1 228 write(std_out,*)" gvec_kss ",ig,"/",ng,gvec_kss(:,ig),test_gvec_kss(:,ig) 229 end if 230 end do 231 ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss") 232 ABI_FREE(test_gvec_kss) 233 234 ! Get important dimension from Hdr_wfk 235 ! Check also the consistency btw Hdr_wfk and Dtset. 236 Ep%nsppol=Hdr_wfk%nsppol 237 Ep%nkibz =Hdr_wfk%nkpt 238 239 if (Ep%nbnds>mband) then 240 Ep%nbnds=mband 241 Dtset%nband(:)=mband 242 write(msg,'(4a,i4,a)')ch10,& 243 & ' Number of bands found less then required. ',ch10,& 244 & ' Calculation will proceed with nbnds = ',mband,ch10 245 MSG_WARNING(msg) 246 end if 247 248 call crystal_from_hdr(Cryst,Hdr_wfk,timrev,remove_inv) 249 call crystal_print(Cryst,mode_paral='COLL') 250 251 ! === Create basic data types for the calculation === 252 ! * Kmesh defines the k-point sampling for the wavefunctions. 253 ! * Qmesh defines the q-point sampling for chi0, all possible differences k1-k2 reduced to the IBZ. 254 ! TODO Kmesh%bz should be [-half,half[ but this modification will be painful! 255 256 call kmesh_init(Kmesh,Cryst,Ep%nkibz,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.) 257 !call kmesh_init(Kmesh,Cryst,Ep%nkibz,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.) 258 ! Some required information are not filled up inside kmesh_init 259 ! So doing it here, even though it is not clean 260 Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:) 261 Kmesh%nshift =Dtset%nshiftk 262 ABI_ALLOCATE(Kmesh%shift,(3,Kmesh%nshift)) 263 Kmesh%shift(:,:) =Dtset%shiftk(:,1:Dtset%nshiftk) 264 265 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 266 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",ab_out, 0, "COLL") 267 268 ! === Find Q-mesh, and do setup for long wavelength limit === 269 ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 270 ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed) 271 call find_qmesh(Qmesh,Cryst,Kmesh) 272 273 call kmesh_print(Qmesh,"Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL") 274 call kmesh_print(Qmesh,"Q-mesh for the screening function",ab_out ,0 ,"COLL") 275 276 do iqbz=1,Qmesh%nbz 277 call get_BZ_item(Qmesh,iqbz,qbz,iq_ibz,isym,itim) 278 sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz)) 279 if (ANY(ABS(qbz-sq )>1.0d-4)) then 280 write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')& 281 & ' qpoint ',qbz,' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,& 282 & ' through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,& 283 & ' however a non zero umklapp G_o vector is required and this is not yet allowed' 284 MSG_ERROR(msg) 285 end if 286 end do 287 288 ! Write the list of qpoints for the screening in netcdf format and exit. 289 ! This file is used by abipy to generate multiple input files. 290 if (Dtset%nqptdm == -1) then 291 if (my_rank==master) then 292 #ifdef HAVE_NETCDF 293 ncerr = nctk_write_ibz(strcat(dtfil%filnam_ds(4), "_qptdms.nc"), qmesh%ibz, qmesh%wt) 294 NCF_CHECK(ncerr) 295 #endif 296 end if 297 MSG_ERROR_NODUMP("Aborting now") 298 end if 299 300 if (Dtset%gw_nqlwl==0) then 301 Ep%nqlwl=1 302 ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl)) 303 Ep%qlwl(:,1)=GW_Q0_DEFAULT ! Use default shift 0.000010, 0.000020, 0.000030 304 else 305 Ep%nqlwl=Dtset%gw_nqlwl 306 ABI_MALLOC(Ep%qlwl,(3,Ep%nqlwl)) 307 Ep%qlwl(:,:)=Dtset%gw_qlwl(:,1:Ep%nqlwl) 308 ABI_CHECK(Ep%nqlwl==1,"nqlwl/=1 not coded") 309 end if 310 !write(std_out,*)" Using qlwl = ",Ep%qlwl 311 312 ! Find optimal value for G-sphere enlargment due to oscillator matrix elements 313 call get_ng0sh(Kmesh%nbz,Kmesh%bz,Qmesh%nibz,Qmesh%ibz,Kmesh%nbz,Kmesh%bz,GW_TOLQ0,ng0sh_opt) 314 call wrtout(std_out,sjoin(' Optimal value for ng0sh:',ltoa(ng0sh_opt)),"COLL") 315 316 Ep%mG0(:)=ng0sh_opt(:) !Ep%mG0(:) = [3, 3, 3] 317 318 ! === In case of symmetrization, find the little group of the q"s === 319 ! * For the long-wavelength limit we consider a small but finite q. However the oscillators are 320 ! evaluated setting q==0. Thus it is possible to take advantage of symmetries also when q --> 0. 321 ! * Here we calculate the enlargement of the G-sphere, npwepG0, needed to account for umklapps. 322 ! TODO Switch on use_umklp, write all this stuff to ab_out 323 324 Ep%npwepG0=Ep%npwe 325 ABI_DT_MALLOC(Ltg_q,(Qmesh%nibz)) 326 327 do iq=1,Qmesh%nibz 328 qtmp(:)=Qmesh%ibz(:,iq); if (normv(qtmp,gmet,'G')<GW_TOLQ0) qtmp(:)=zero; use_umklp=0 329 call littlegroup_init(qtmp,Kmesh,Cryst,use_umklp,Ltg_q(iq),Ep%npwe,gvec=gvec_kss) 330 end do 331 332 ecutepspG0 = Dtset%ecuteps 333 ABI_CHECK(ecutepspG0 > zero, "ecuteps must be > 0") 334 if (Ep%symchi/=0) then 335 ecutepspG0=MAXVAL(Ltg_q(:)%max_kin_gmG0)+tol6; npwepG0=0; nshepspG0=0 336 write(std_out,*)" Due to umklapp processes : ecutepspg0= ",ecutepspG0 337 call setshells(ecutepspG0,npwepG0,nshepspG0,Cryst%nsym,gmet,gprimd,Cryst%symrel,'eps_pG0',Cryst%ucvol) 338 Ep%npwepG0=npwepG0 339 end if 340 341 if (Ep%npwepG0>Ep%npwvec) then 342 write(msg,'(3a,i5,a,i5)')& 343 & ' npwepG0 > npwvec, decrease npweps or increase npwwfn. ',ch10,& 344 & ' npwepG0 = ',Ep%npwepG0,' npwvec = ',Ep%npwvec 345 MSG_ERROR(msg) 346 end if 347 348 ! === Create structure describing the G-sphere used for chi0/espilon and Wfns === 349 ! * The cutoff is >= ecuteps to allow for umklapp 350 #if 0 351 call gsph_init(Gsph_wfn,Cryst,0,ecut=Dtset%ecutwfn) 352 Dtset%npwwfn = Gsph_wfn%ng 353 Ep%npwwfn = Gsph_wfn%ng 354 ierr = 0 355 do ig=1,MIN(Gsph_wfn%ng, ng_kss) 356 if ( ANY(Gsph_wfn%gvec(:,ig) /= gvec_kss(:,ig)) ) then 357 write(std_out,*)ig, Gsph_wfn%gvec(:,ig), gvec_kss(:,ig) 358 end if 359 end do 360 ABI_CHECK(ierr==0,"Wrong gvec_wfn") 361 #else 362 call gsph_init(Gsph_wfn,Cryst,Ep%npwvec,gvec=gvec_kss) 363 #endif 364 365 #if 0 366 call gsph_init(Gsph_epsG0,Cryst,0,ecut=ecutepspG0) 367 Ep%npwepG0 = Gsph_epsG0%ng 368 ierr = 0 369 do ig=1,MIN(Gsph_epsG0%ng, ng_kss) 370 if ( ANY(Gsph_epsG0%gvec(:,ig) /= gvec_kss(:,ig)) ) then 371 write(std_out,*)ig, Gsph_epsG0%gvec(:,ig), gvec_kss(:,ig) 372 end if 373 end do 374 ABI_CHECK(ierr==0,"Wrong gvec_epsG0") 375 #else 376 call gsph_init(Gsph_epsG0,Cryst,Ep%npwepG0,gvec=gvec_kss) 377 #endif 378 ! 379 ! ======================================================================= 380 ! ==== Setup of the FFT mesh used for the oscillator matrix elements ==== 381 ! ======================================================================= 382 ! * ngfft_gw(7:18) is the same as Dtset%ngfft(7:18), initialized before entering setup_screening. 383 ! Here we just redefine ngfft_gw(1:6) according to the following options: 384 ! 385 ! method==0 ==> FFT grid read from __fft.in__ (only for debugging purpose) 386 ! method==1 ==> normal FFT grid 387 ! method==2 ==> slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 388 ! method==3 ==> doubled FFT grid, to treat exactly the convolution defining the density, 389 ! Useful in sigma if ppmodel=[2,3,4] since rho(G-Gp) or to calculate matrix elements of v_Hxc. 390 ! 391 ! enforce_sym==1 ==> enforce a direct space FFT mesh compatible with all symmetries operation 392 ! enforce_sym==0 ==> Find the smallest FFT grid compatibile with the library, do not care about symmetries 393 ! 394 ngfft_gw(1:18)=Dtset%ngfft(1:18); method=2 395 if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0 396 if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1 397 if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2 398 if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3 399 enforce_sym=MOD(Dtset%fftgw,10) 400 401 ! Use npwepG0 to account for umklapps. 402 call setmesh(gmet,gvec_kss,ngfft_gw,Ep%npwvec,Ep%npwepG0,Ep%npwwfn,nfftgw_tot,method,Ep%mG0,Cryst,enforce_sym) 403 !call new_setmesh(Cryst,ecut_osc,ecutwfn,nkpt,kpoints,method,Ep%mG0,enforce_sym,ngfft_gw,nfftgw_tot) 404 405 ABI_FREE(gvec_kss) 406 407 ! FIXME this wont work if nqptdm/=0 408 call vcoul_init(Vcp,Gsph_epsG0,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecuteps,Ep%npwe,Ep%nqlwl,& 409 & Ep%qlwl,ngfftf,comm) 410 411 #if 0 412 ! Using the random q for the optical limit is one of the reasons 413 ! why sigma breaks the initial energy degeneracies. 414 Vcp%i_sz=zero 415 Vcp%vc_sqrt(1,:)=czero 416 Vcp%vcqlwl_sqrt(1,:)=czero 417 #endif 418 419 ! Value of scissor energy 420 Ep%mbpt_sciss=Dtset%mbpt_sciss 421 422 ! Define the frequency mesh for epsilon according to the method used. 423 Ep%nomegaei=1 424 Ep%nomegaer=1; if (is_static) Ep%nomegaer=0 425 Ep%nomegaec=0 426 Ep%omegaermax=zero 427 428 ! For ppmodels 2,3,4, only omega=0 is needed. 429 if (Ep%plasmon_pole_model.and.Dtset%nfreqre==1.and.Dtset%nfreqim==0) then 430 Ep%nomegaer=1; Ep%nomegaei=0 431 write(msg,'(7a)')ch10,& 432 & ' The inverse dielectric matrix will be calculated on zero frequency only',ch10,& 433 & ' please note that the calculated epsilon^-1 cannot be used ',ch10,& 434 & ' to calculate QP corrections using plasmonpole model 1',ch10 435 call wrtout(std_out,msg,'COLL') 436 call wrtout(ab_out,msg,'COLL') 437 end if 438 439 ! Max number of omega along the imaginary axis 440 if (Ep%analytic_continuation.or.Ep%contour_deformation) then 441 Ep%nomegaei=Dtset%nfreqim 442 if (Dtset%gw_frqim_inzgrid==1) then 443 MSG_WARNING('iomega = z/1-z transfom grid will be used for imaginary frequency grid') 444 end if 445 if (Dtset%cd_customnimfrqs/=0) then 446 MSG_WARNING('Custom imaginary grid specified. Assuming experienced user.') 447 Ep%nomegaei=Dtset%cd_customnimfrqs 448 end if 449 if (Ep%nomegaei==-1) then 450 Ep%nomegaei=NOMEGAGAUSS 451 MSG_WARNING(sjoin('Number of imaginary frequencies set to default= ',itoa(NOMEGAGAUSS))) 452 end if 453 if (Ep%nomegaei==0) then 454 MSG_WARNING(' nfreqim = 0! Assuming experienced user merging several frequency calculations.') 455 end if 456 end if 457 458 ! Range and total number of real frequencies. 459 Ep%omegaermin = zero 460 if (Ep%contour_deformation) then 461 Ep%nomegaer=Dtset%nfreqre; Ep%omegaermin=Dtset%freqremin; Ep%omegaermax=Dtset%freqremax 462 if (Dtset%gw_frqre_tangrid==1) then 463 Ep%omegaermax=Dtset%cd_max_freq 464 MSG_WARNING('Tangent transfom grid will be used for real frequency grid') 465 end if 466 if (Dtset%gw_frqre_tangrid==1) then 467 MSG_WARNING('Tangent transfom grid will be used for real frequency grid') 468 end if 469 if (Ep%nomegaer==-1) then 470 Ep%nomegaer=NOMEGAREAL 471 MSG_WARNING(sjoin('Number of real frequencies set to default= ',itoa(NOMEGAREAL))) 472 end if 473 if (Ep%nomegaer==0) then 474 MSG_WARNING('nfreqre = 0 ! Assuming experienced user.') 475 end if 476 if (ABS(Ep%omegaermin)<TOL16) then 477 Ep%omegaermin=zero 478 write(msg,'(a,f8.4)')' Min real frequency set to default [Ha] = ',Ep%omegaermin 479 MSG_WARNING(msg) 480 end if 481 if (Ep%omegaermin>Ep%omegaermax) then 482 MSG_ERROR('freqremin > freqremax !') 483 end if 484 if (Ep%omegaermax<TOL16) then 485 Ep%omegaermax=OMEGAERMAX 486 write(msg,'(a,f8.4)')' Max real frequency set to default [Ha] = ',OMEGAERMAX 487 MSG_WARNING(msg) 488 end if 489 ! Check if a subset of the frequencies is to be used 490 if (Dtset%cd_subset_freq(1)/=0) then 491 istart = Dtset%cd_subset_freq(1) 492 iend = Dtset%cd_subset_freq(2) 493 if (istart>iend.or.istart<0.or.iend<0) then 494 MSG_ERROR(' check indices of cd_subset_freq!') 495 end if 496 write(msg,'(2(a,i0))')' Using cd_subset_freq to only do freq. from ',istart,' to ',iend 497 MSG_WARNING(msg) 498 ! Reset the numbers 499 if (Dtset%gw_frqre_tangrid/=1) then ! Normal equidistant grid 500 Ep%nomegaer = iend-istart+1 501 domegareal=(Ep%omegaermax-Ep%omegaermin)/(Ep%nomegaer-1) 502 Ep%omegaermin = Ep%omegaermin+(istart-1)*domegareal 503 Ep%omegaermax = Ep%omegaermin+(iend-1)*domegareal 504 else 505 Ep%nomegaer = iend-istart+1 506 end if 507 end if 508 end if 509 510 ! Check full grid calculations 511 if (Dtset%cd_full_grid/=0) then 512 MSG_WARNING("FULL GRID IN COMPLEX PLANE CALCULATED. YOU MIGHT NOT BE ABLE TO USE SCREENING FILES!") 513 if (Dtset%cd_subset_freq(1)/=0) then 514 MSG_ERROR('cd_subset_freq cannot be used with cd_full_grid!') 515 end if 516 Ep%nomegaec = Ep%nomegaei*(Ep%nomegaer-1) 517 end if 518 519 Ep%nomega=Ep%nomegaer+Ep%nomegaei+Ep%nomegaec ! Total number of frequencies. 520 521 ! ==== Setup of the spectral method ==== 522 Ep%spmeth =Dtset%spmeth; Ep%nomegasf=Dtset%nomegasf; Ep%spsmear =Dtset%spbroad 523 524 if (Ep%spmeth/=0) then 525 write(msg,'(2a,i3,2a,i8)')ch10,& 526 & ' setup_screening: using spectral method: ',Ep%spmeth,ch10,& 527 & ' Number of frequencies for imaginary part: ',Ep%nomegasf 528 call wrtout(std_out,msg,'COLL') 529 if (Ep%spmeth==2) then 530 write(msg,'(a,f8.5,a)')' Gaussian broadening = ',Ep%spsmear*Ha_eV,' [eV]' 531 call wrtout(std_out,msg,'COLL') 532 end if 533 end if 534 535 Ep%nI=1; Ep%nJ=1 536 if (Dtset%nspinor==2) then 537 !if (Dtset%usepaw==1.and.Dtset%pawspnorb>0) then 538 ! Ep%nI=1; Ep%nJ=4 539 !end if 540 ! For spin-spin interaction 541 ! Ep%nI=4; Ep%nJ=4 542 ABI_CHECK(Ep%npwepG0 == Ep%npwe, "npwepG0 must be == npwe if spinor==2") 543 !ABI_CHECK(Ep%symchi == 0, "symchi/=0 and nspinor=2 not available") 544 end if 545 546 ! === Enable the calculations of chi0 on user-specified q-points === 547 Ep%nqibz=Qmesh%nibz 548 ABI_MALLOC(Ep%qibz,(3,Ep%nqibz)) 549 Ep%qibz(:,:)=Qmesh%ibz(:,:) 550 551 Ep%nqcalc=Ep%nqibz 552 if (Dtset%nqptdm>0) Ep%nqcalc=Dtset%nqptdm 553 554 ABI_MALLOC(Ep%qcalc,(3,Ep%nqcalc)) 555 if (Ep%nqcalc/=Ep%nqibz) then 556 write(msg,'(6a)')ch10,& 557 & ' Dielectric matrix will be calculated only for some ',ch10,& 558 & ' selected q points provided by the user through the input variables ',ch10,& 559 & ' nqptdm and qptdm' 560 call wrtout(std_out,msg,'COLL') 561 call wrtout(ab_out,msg,'COLL') 562 ltest= Ep%nqcalc <= Qmesh%nibz 563 ABI_CHECK(ltest, 'nqptdm should not exceed the number of q points in the IBZ') 564 Ep%qcalc(:,:)=Dtset%qptdm(:,1:Ep%nqcalc) 565 ! Check whether the q-points provided are correct. 566 do iq=1,Ep%nqcalc 567 found=.FALSE. 568 do iqp=1,Qmesh%nibz 569 qtmp(:)=Ep%qcalc(:,iq)-Qmesh%ibz(:,iqp) 570 found=(normv(qtmp,gmet,'G')<GW_TOLQ) 571 if (found) EXIT 572 end do 573 ABI_CHECK(found, 'One or more points specified by Dtset%qptdm do not satisfy q=k1-k2') 574 end do 575 else 576 Ep%qcalc(:,:)=Ep%qibz(:,:) 577 end if 578 579 ! To write the SCR header correctly, with heads and wings, we have 580 ! to make sure that q==0, if present, is the first q-point in the list. 581 !has_q0=(ANY(normv(Ep%qcalc(:,:),gmet,'G')<GW_TOLQ0)) !commented to avoid problems with sunstudio12 582 has_q0=.FALSE. 583 do iq=1,Ep%nqcalc 584 if (normv(Ep%qcalc(:,iq),gmet,'G')<GW_TOLQ0) then 585 has_q0=.TRUE.; EXIT 586 end if 587 end do 588 589 if (has_q0.and.normv(Ep%qcalc(:,1),gmet,'G')>=GW_TOLQ0) then 590 write(msg,'(5a)')& 591 & 'The list of q-points to be calculated contains the Gamma point, ',ch10,& 592 & 'however Gamma is not the first point in the list. ' ,ch10,& 593 & 'Please, change your input file accordingly. ' 594 MSG_ERROR(msg) 595 end if 596 597 ! === Initialize the band structure datatype === 598 ! * Copy KSS energies and occupations up to Ep%nbnds==Dtset%nband(:) 599 ! TODO Recheck symmorphy and inversion 600 bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)) 601 602 ABI_MALLOC(doccde,(bantot)) 603 ABI_MALLOC(eigen,(bantot)) 604 ABI_MALLOC(occfact,(bantot)) 605 doccde(:)=zero; eigen(:)=zero; occfact(:)=zero 606 607 jj=0; ibtot=0 608 do isppol=1,Dtset%nsppol 609 do ikibz=1,Dtset%nkpt 610 do ib=1,Hdr_wfk%nband(ikibz+Dtset%nkpt*(isppol-1)) 611 ibtot=ibtot+1 612 if (ib<=Ep%nbnds) then 613 jj=jj+1 614 occfact(jj)=Hdr_wfk%occ(ibtot) 615 eigen (jj)=energies_p(ib,ikibz,isppol) 616 end if 617 end do 618 end do 619 end do 620 ABI_FREE(energies_p) 621 622 ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of 623 ! the symmetry operations in the old GW code (symmorphy and inversion) 624 ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6)) 625 ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt') 626 627 ABI_MALLOC(npwarr,(Hdr_wfk%nkpt)) 628 npwarr(:)=Ep%npwwfn 629 630 call ebands_init(bantot,KS_BSt,Dtset%nelect,doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,& 631 & Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,& 632 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 633 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 634 635 ! TODO modify outkss in order to calculate the eigenvalues also if NSCF calculation. 636 ! this fails simply because in case of NSCF occ are zero 637 !ltest=(ALL(ABS(occfact-KS_BSt%occ)<1.d-2)) 638 !call assert(ltest,'difference in occfact') 639 !write(std_out,*)MAXVAL(ABS(occfact(:)-KS_BSt%occ(:))) 640 641 !TODO call ebands_update_occ here 642 !$call ebands_update_occ(KS_BSt,spinmagntarget,stmbias,Dtset%prtvol) 643 644 ABI_FREE(doccde) 645 ABI_FREE(eigen) 646 ABI_FREE(npwarr) 647 648 ! Initialize abinit header for the screening part 649 call hdr_init(KS_BSt,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl) 650 651 ! Get Pawrhoij from the header. 652 ABI_DT_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw)) 653 if (Dtset%usepaw==1) then 654 call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 655 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 656 end if 657 call hdr_update(Hdr_out,bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 658 659 ABI_FREE(occfact) 660 call pawrhoij_free(Pawrhoij) 661 ABI_DT_FREE(Pawrhoij) 662 663 ! ==== Setup of extrapolar technique ==== 664 Ep%gwcomp = Dtset%gwcomp; Ep%gwencomp = Dtset%gwencomp 665 if (Ep%gwcomp == 1) then 666 write(msg,'(a,f8.2,a)')' Using the completeness correction with gwencomp ',Ep%gwencomp*Ha_eV,' [eV] ' 667 call wrtout(std_out,msg,'COLL') 668 end if 669 670 ! Final compatibility tests 671 if (ANY(KS_BSt%istwfk /= 1)) then 672 MSG_WARNING('istwfk/=1 is still under development') 673 end if 674 675 ltest = (KS_BSt%mband == Ep%nbnds .and. ALL(KS_BSt%nband == Ep%nbnds)) 676 ABI_CHECK(ltest,'BUG in definition of KS_BSt%nband') 677 678 if (Ep%gwcomp==1 .and. Ep%spmeth>0) then 679 MSG_ERROR("Hilbert transform and extrapolar method are not compatible") 680 end if 681 682 DBG_EXIT('COLL') 683 684 end subroutine setup_screening