TABLE OF CONTENTS
ABINIT/setup_bse [ Functions ]
NAME
setup_bse
FUNCTION
This routine performs the initialization of basic objects and quantities used for Bethe-Salpeter calculations. In particular the excparam data type that defines the parameters of the calculation is completely initialized starting from the content of Dtset and the parameters read from the external WFK and SCR (SUSC) file.
COPYRIGHT
Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida) Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
codvsn=Code version ngfft_gw(18)=Information about 3D FFT for density and potentials, see ~abinit/doc/variables/vargs.htm#ngfft acell(3)=Length scales of primitive translations (bohr) rprim(3,3)=Dimensionless real space primitive translations. Dtset<dataset_type>=All input variables for this dataset. Some of them might be redefined here TODO Dtfil=filenames and unit numbers used in abinit. Psps <pseudopotential_type>=variables related to pseudopotentials Pawtab(Psps%ntypat*Dtset%usepaw)<pawtab_type>=PAW tabulated starting data
OUTPUT
Cryst<crystal_t>=Info on the crystalline Structure. Kmesh<kmesh_t>=Structure defining the k-sampling for the wavefunctions. Qmesh<kmesh_t>=Structure defining the q-sampling for the symmetrized inverse dielectric matrix. Gsph_x<gsphere_t=Data type gathering info on the G-sphere for wave functions and e^{-1}, KS_BSt<ebands_t>=The KS band structure (energies, occupancies, k-weights...) Vcp<vcoul_t>=Structure gathering information on the Coulomb interaction in reciprocal space, including a possible cutoff in real space. ngfft_osc(18)=Contain all needed information about the 3D FFT for the oscillator matrix elements. See ~abinit/doc/variables/vargs.htm#ngfft Bsp<excparam>=Basic parameters defining the Bethe-Salpeter run. Completely initialed in output. Hdr_wfk<Hdr_type>=The header of the WFK file. Hdr_bse<Hdr_type>=Local header initialized from the parameters used for the Bethe-Salpeter calculation. BS_files<excfiles>=Files used in the calculation. w_file=File name used to construct W. Set to ABI_NOFILE if no external file is used.
PARENTS
bethe_salpeter
CHILDREN
apply_scissor,bsp_calctype2str,crystal_from_hdr,crystal_print ebands_copy,ebands_init,ebands_print,ebands_report_gap ebands_update_occ,eprenorms_bcast,eprenorms_from_epnc,find_qmesh get_bz_item,get_ng0sh,gsph_extend,gsph_init,hdr_init,hdr_update hdr_vs_dtset,hscr_bcast,hscr_free,hscr_from_file,hscr_print init_transitions,kmesh_init,kmesh_print,make_mesh,matrginv,metric mkrdim,pawrhoij_alloc,pawrhoij_copy,pawrhoij_free,print_bs_files print_bs_parameters,print_gsphere,print_ngfft,rdgw,setmesh,vcoul_init wfk_read_eigenvalues,wrtout,xmpi_bcast,xmpi_max,xmpi_split_work
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 subroutine setup_bse(codvsn,acell,rprim,ngfftf,ngfft_osc,Dtset,Dtfil,BS_files,Psps,Pawtab,BSp,& 69 & Cryst,Kmesh,Qmesh,KS_BSt,QP_bst,Hdr_wfk,Gsph_x,Gsph_c,Vcp,Hdr_bse,w_fname,Epren,comm,Wvl) 70 71 use defs_basis 72 use defs_datatypes 73 use defs_abitypes 74 use defs_wvltypes 75 use m_bs_defs 76 use m_profiling_abi 77 use m_errors 78 use m_xmpi 79 use m_nctk 80 use m_hdr 81 82 use m_gwdefs, only : GW_Q0_DEFAULT 83 use m_fstrings, only : toupper, sjoin 84 use m_io_tools, only : file_exists, open_file 85 use m_geometry, only : normv, mkrdim, metric 86 use m_abilasi, only : matrginv 87 use m_crystal, only : crystal_print, idx_spatial_inversion, crystal_t 88 use m_crystal_io, only : crystal_from_hdr 89 use m_bz_mesh, only : kmesh_t, kmesh_init, get_ng0sh, kmesh_print, get_BZ_item, find_qmesh, make_mesh 90 use m_ebands, only : ebands_init, ebands_print, ebands_copy, ebands_free, & 91 & ebands_update_occ, get_valence_idx, apply_scissor, ebands_report_gap 92 use m_eprenorms, only : eprenorms_t, eprenorms_from_epnc, eprenorms_bcast 93 use m_vcoul, only : vcoul_t, vcoul_init 94 use m_fftcore, only : print_ngfft 95 use m_fft_mesh, only : setmesh 96 use m_gsphere, only : gsphere_t, gsph_init, print_gsphere, merge_and_sort_kg, gsph_extend 97 use m_io_screening, only : hscr_t, hscr_free, hscr_io, hscr_bcast, hscr_from_file, hscr_print 98 use m_pawtab, only : pawtab_type 99 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free 100 use m_qparticles, only : rdgw 101 use m_screen, only : MDL_BECHSTEDT 102 use m_wfk, only : wfk_read_eigenvalues 103 104 !This section has been created automatically by the script Abilint (TD). 105 !Do not modify the following lines by hand. 106 #undef ABI_FUNC 107 #define ABI_FUNC 'setup_bse' 108 use interfaces_14_hidewrite 109 use interfaces_56_io_mpi 110 !End of the abilint section 111 112 implicit none 113 114 !Arguments ------------------------------------ 115 !scalars 116 integer,intent(in) :: comm 117 character(len=6),intent(in) :: codvsn 118 character(len=fnlen),intent(out) :: w_fname 119 type(dataset_type),intent(inout) :: Dtset 120 type(datafiles_type),intent(in) :: Dtfil 121 type(pseudopotential_type),intent(in) :: Psps 122 type(excparam),intent(inout) :: Bsp 123 type(hdr_type),intent(out) :: Hdr_wfk,Hdr_bse 124 type(crystal_t),intent(out) :: Cryst 125 type(kmesh_t),intent(out) :: Kmesh,Qmesh 126 type(gsphere_t),intent(out) :: Gsph_x,Gsph_c 127 type(ebands_t),intent(out) :: KS_BSt,QP_Bst 128 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw) 129 type(vcoul_t),intent(out) :: Vcp 130 type(excfiles),intent(out) :: BS_files 131 type(wvl_internal_type), intent(in) :: Wvl 132 type(eprenorms_t),intent(out) :: Epren 133 !arrays 134 integer,intent(in) :: ngfftf(18) 135 integer,intent(out) :: ngfft_osc(18) 136 real(dp),intent(in) :: acell(3),rprim(3,3) 137 138 !Local variables ------------------------------ 139 !scalars 140 integer,parameter :: pertcase0=0,master=0 141 integer(i8b) :: work_size,tot_nreh,neh_per_proc,il 142 integer :: bantot,enforce_sym,ib,ibtot,ik_ibz,isppol,jj,method,iat,ount !ii, 143 integer :: mband,io,nfftot_osc,spin,hexc_size,nqlwl,iq 144 integer :: timrev,iq_bz,isym,iq_ibz,itim 145 integer :: my_rank,nprocs,fform,npwe_file,ierr,my_k1, my_k2,my_nbks 146 integer :: first_dig,second_dig,it 147 real(dp) :: ucvol,qnorm 148 real(dp):: eff,mempercpu_mb,wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem 149 logical,parameter :: remove_inv=.FALSE. 150 logical :: ltest,occ_from_dtset 151 character(len=500) :: msg 152 character(len=fnlen) :: gw_fname,test_file,wfk_fname 153 character(len=fnlen) :: ep_nc_fname 154 type(hscr_t) :: Hscr 155 !arrays 156 integer :: ng0sh_opt(3),val_idx(Dtset%nsppol) 157 integer,allocatable :: npwarr(:),val_indeces(:,:),nlmn_atm(:) 158 real(dp) :: qpt_bz(3),minmax_tene(2) 159 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3) 160 real(dp) :: qred2cart(3,3),qcart2red(3,3) 161 real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:) 162 real(dp),allocatable :: igwene(:,:,:) 163 real(dp),pointer :: energies_p(:,:,:) 164 complex(dpc),allocatable :: gw_energy(:,:,:) 165 type(Pawrhoij_type),allocatable :: Pawrhoij(:) 166 !Interp@BSE 167 !integer :: mode 168 !integer :: kmult(3) 169 !integer :: unt 170 !integer :: rl_nb 171 !logical :: interp_params_exists, prepare, sum_overlaps 172 !namelist /interp_params/ mode,kmult,prepare,rl_nb,sum_overlaps 173 !character(len=fnlen) :: tmp_fname 174 175 !************************************************************************ 176 177 DBG_ENTER("COLL") 178 179 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 180 181 ! === Check for calculations that are not implemented === 182 ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1)) 183 ABI_CHECK(ltest,'Dtset%nband must be constant') 184 ABI_CHECK(Dtset%nspinor==1,"nspinor==2 not coded") 185 186 ! === Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume === 187 call mkrdim(acell,rprim,rprimd) 188 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 189 190 ! Read energies and header from the WFK file. 191 wfk_fname = dtfil%fnamewffk 192 if (.not. file_exists(wfk_fname)) then 193 wfk_fname = nctk_ncify(wfk_fname) 194 MSG_COMMENT(sjoin("File not found. Will try netcdf file: ", wfk_fname)) 195 end if 196 197 call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm) 198 mband = MAXVAL(Hdr_wfk%nband) 199 200 call hdr_vs_dtset(Hdr_wfk,Dtset) 201 202 ! === Create crystal_t data type === 203 !remove_inv= .FALSE. !(nsym_kss/=Hdr_wfk%nsym) 204 timrev= 2 ! This information is not reported in the header 205 ! 1 => do not use time-reversal symmetry 206 ! 2 => take advantage of time-reversal symmetry 207 208 call crystal_from_hdr(Cryst,Hdr_wfk,timrev,remove_inv) 209 call crystal_print(Cryst) 210 ! 211 ! Setup of the k-point list and symmetry tables in the BZ ----------------------------------- 212 if (Dtset%chksymbreak==0) then 213 MSG_WARNING("Calling make_mesh") 214 call make_mesh(Kmesh,Cryst,Dtset%kptopt,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk,break_symmetry=.TRUE.) 215 ! TODO 216 !Check if kibz from KSS file corresponds to the one returned by make_mesh. 217 else 218 call kmesh_init(Kmesh,Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt) 219 end if 220 BSp%nkibz = Kmesh%nibz !We might allow for a smaller number of points.... 221 222 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 223 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",ab_out, 0, "COLL") 224 225 nqlwl = 0 226 w_fname = ABI_NOFILE 227 if (Dtset%getscr/=0.or.Dtset%irdscr/=0) then 228 w_fname=Dtfil%fnameabi_scr 229 else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then 230 w_fname=Dtfil%fnameabi_sus 231 MSG_ERROR("(get|ird)suscep not implemented") 232 end if 233 234 if (w_fname /= ABI_NOFILE) then ! Read dimensions from the external file 235 236 if (.not. file_exists(w_fname)) then 237 w_fname = nctk_ncify(w_fname) 238 MSG_COMMENT(sjoin("File not found. Will try netcdf file: ", w_fname)) 239 end if 240 241 if (my_rank==master) then 242 ! Master reads npw and nqlwl from SCR file. 243 call wrtout(std_out,sjoin('Testing file: ', w_fname),"COLL") 244 245 call hscr_from_file(hscr,w_fname,fform,xmpi_comm_self) 246 ! Echo the header. 247 if (Dtset%prtvol>0) call hscr_print(Hscr) 248 249 npwe_file = Hscr%npwe ! Have to change %npweps if it was larger than dim on disk. 250 nqlwl = Hscr%nqlwl 251 252 if (Dtset%npweps>npwe_file) then 253 write(msg,'(2(a,i0),2a,i0)')& 254 & "The number of G-vectors stored on file (",npwe_file,") is smaller than Dtset%npweps = ",Dtset%npweps,ch10,& 255 & "Calculation will proceed with the maximum available set, npwe_file = ",npwe_file 256 MSG_WARNING(msg) 257 Dtset%npweps = npwe_file 258 else if (Dtset%npweps<npwe_file) then 259 write(msg,'(2(a,i0),2a,i0)')& 260 & "The number of G-vectors stored on file (",npwe_file,") is larger than Dtset%npweps = ",Dtset%npweps,ch10,& 261 & "Calculation will proceed with Dtset%npweps = ",Dtset%npweps 262 MSG_COMMENT(msg) 263 end if 264 end if 265 266 call hscr_bcast(Hscr,master,my_rank,comm) 267 call xmpi_bcast(Dtset%npweps,master,comm,ierr) 268 call xmpi_bcast(nqlwl,master,comm,ierr) 269 270 if (nqlwl>0) then 271 ABI_MALLOC(qlwl,(3,nqlwl)) 272 qlwl = Hscr%qlwl 273 end if 274 ! 275 ! Init Qmesh from the SCR file. 276 call kmesh_init(Qmesh,Cryst,Hscr%nqibz,Hscr%qibz,Dtset%kptopt) 277 278 ! The G-sphere for W and Sigma_c is initialized from the gvectors found in the SCR file. 279 call gsph_init(Gsph_c,Cryst,Dtset%npweps,gvec=Hscr%gvec) 280 281 call hscr_free(Hscr) 282 else 283 ! Init Qmesh from the K-mesh reported in the WFK file. 284 call find_qmesh(Qmesh,Cryst,Kmesh) 285 286 ! The G-sphere for W and Sigma_c is initialized from ecuteps. 287 call gsph_init(Gsph_c,Cryst,0,ecut=Dtset%ecuteps) 288 Dtset%npweps = Gsph_c%ng 289 end if 290 291 BSp%npweps = Dtset%npweps 292 BSp%ecuteps = Dtset%ecuteps 293 294 if (nqlwl==0) then 295 nqlwl=1 296 ABI_MALLOC(qlwl,(3,nqlwl)) 297 qlwl(:,nqlwl)= GW_Q0_DEFAULT 298 write(msg,'(3a,i2,a,3f9.6)')& 299 & "The Header of the screening file does not contain the list of q-point for the optical limit ",ch10,& 300 & "Using nqlwl= ",nqlwl," and qlwl = ",qlwl(:,1) 301 MSG_COMMENT(msg) 302 end if 303 write(std_out,*)"nqlwl and qlwl for Coulomb singularity and e^-1",nqlwl,qlwl 304 305 ! === Setup of q-mesh in the whole BZ === 306 ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 307 ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed) 308 ! 309 call kmesh_print(Qmesh,"Q-mesh for the screening function",std_out,Dtset%prtvol,"COLL") 310 call kmesh_print(Qmesh,"Q-mesh for the screening function",ab_out ,0 ,"COLL") 311 312 do iq_bz=1,Qmesh%nbz 313 call get_BZ_item(Qmesh,iq_bz,qpt_bz,iq_ibz,isym,itim) 314 sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz)) 315 if (ANY(ABS(Qmesh%bz(:,iq_bz)-sq )>1.0d-4)) then 316 write(std_out,*) sq,Qmesh%bz(:,iq_bz) 317 write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')& 318 & 'qpoint ',Qmesh%bz(:,iq_bz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,& 319 & 'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,& 320 & 'however a non zero umklapp G_o vector is required and this is not yet allowed' 321 MSG_ERROR(msg) 322 end if 323 end do 324 325 BSp%algorithm = Dtset%bs_algorithm 326 BSp%nstates = Dtset%bs_nstates 327 Bsp%nsppol = Dtset%nsppol 328 Bsp%hayd_term = Dtset%bs_hayd_term 329 ! 330 ! Define the algorithm for solving the BSE. 331 if (BSp%algorithm == BSE_ALGO_HAYDOCK) then 332 BSp%niter = Dtset%bs_haydock_niter 333 BSp%haydock_tol = Dtset%bs_haydock_tol 334 335 else if (BSp%algorithm == BSE_ALGO_CG) then 336 ! FIXME For the time being use an hardcoded value. 337 ! TODO change name in Dtset% 338 BSp%niter = Dtset%nstep !100 339 BSp%cg_tolwfr = Dtset%tolwfr 340 BSp%nline = Dtset%nline 341 BSp%nbdbuf = Dtset%nbdbuf 342 BSp%nstates = Dtset%bs_nstates 343 MSG_WARNING("Check CG setup") 344 else 345 !BSp%niter = 0 346 !BSp%tol_iter = HUGE(one) 347 end if 348 ! 349 ! Shall we include Local field effects? 350 SELECT CASE (Dtset%bs_exchange_term) 351 CASE (0,1) 352 BSp%exchange_term = Dtset%bs_exchange_term 353 CASE DEFAULT 354 write(msg,'(a,i0)')" Wrong bs_exchange_term: ",Dtset%bs_exchange_term 355 MSG_ERROR(msg) 356 END SELECT 357 ! 358 ! Treatment of the off-diagonal coupling block. 359 SELECT CASE (Dtset%bs_coupling) 360 CASE (0) 361 BSp%use_coupling = 0 362 msg = 'RESONANT ONLY CALCULATION' 363 CASE (1) 364 BSp%use_coupling = 1 365 msg = ' RESONANT+COUPLING CALCULATION ' 366 CASE DEFAULT 367 write(msg,'(a,i0)')" Wrong bs_coupling: ",Dtset%bs_coupling 368 MSG_ERROR(msg) 369 END SELECT 370 call wrtout(std_out,msg,"COLL") 371 372 BSp%use_diagonal_Wgg = .FALSE. 373 Bsp%use_coulomb_term = .TRUE. 374 BSp%eps_inf=zero 375 Bsp%mdlf_type=0 376 377 first_dig =MOD(Dtset%bs_coulomb_term,10) 378 second_dig=Dtset%bs_coulomb_term/10 379 380 Bsp%wtype = second_dig 381 SELECT CASE (second_dig) 382 CASE (BSE_WTYPE_NONE) 383 call wrtout(std_out,"Coulomb term won't be calculated","COLL") 384 Bsp%use_coulomb_term = .FALSE. 385 386 CASE (BSE_WTYPE_FROM_SCR) 387 call wrtout(std_out,"W is read from an external SCR file","COLL") 388 Bsp%use_coulomb_term = .TRUE. 389 390 CASE (BSE_WTYPE_FROM_MDL) 391 call wrtout(std_out,"W is approximated with the model dielectric function","COLL") 392 Bsp%use_coulomb_term = .TRUE. 393 BSp%mdlf_type = MDL_BECHSTEDT 394 BSp%eps_inf = Dtset%mdf_epsinf 395 ABI_CHECK(Bsp%eps_inf > zero, "mdf_epsinf <= 0") 396 397 CASE DEFAULT 398 write(msg,'(a,i0)')" Wrong second digit in bs_coulomb_term: ",Dtset%bs_coulomb_term 399 MSG_ERROR(msg) 400 END SELECT 401 ! 402 ! Diagonal approximation or full matrix? 403 BSp%use_diagonal_Wgg = .TRUE. 404 if (Bsp%wtype /= BSE_WTYPE_NONE) then 405 SELECT CASE (first_dig) 406 CASE (0) 407 call wrtout(std_out,"Using diagonal approximation W_GG","COLL") 408 BSp%use_diagonal_Wgg = .TRUE. 409 CASE (1) 410 call wrtout(std_out,"Using full W_GG' matrix ","COLL") 411 BSp%use_diagonal_Wgg = .FALSE. 412 CASE DEFAULT 413 write(msg,'(a,i0)')" Wrong first digit in bs_coulomb_term: ",Dtset%bs_coulomb_term 414 MSG_ERROR(msg) 415 END SELECT 416 end if 417 418 !TODO move the initialization of the parameters for the interpolation in setup_bse_interp 419 420 BSp%use_interp = .FALSE. 421 BSp%interp_mode = BSE_INTERP_YG 422 BSp%interp_kmult(1:3) = 0 423 BSp%prep_interp = .FALSE. 424 BSp%sum_overlaps = .TRUE. ! Sum over the overlaps 425 426 ! Printing ncham 427 BSp%prt_ncham = .FALSE. 428 429 ! Deactivate Interpolation Technique by default 430 ! if (.FALSE.) then 431 432 ! Reading parameters from the input file 433 BSp%use_interp = (dtset%bs_interp_mode /= 0) 434 BSp%prep_interp = (dtset%bs_interp_prep == 1) 435 436 SELECT CASE (dtset%bs_interp_mode) 437 CASE (0) 438 ! No interpolation, do not print anything ! 439 CASE (1) 440 call wrtout(std_out,"Using interpolation technique with energies and wavefunctions from dense WFK","COLL") 441 CASE (2) 442 call wrtout(std_out,"Interpolation technique with energies and wfn on dense WFK + treatment ABC of divergence","COLL") 443 CASE (3) 444 call wrtout(std_out,"Interpolation technique + divergence ABC along diagonal","COLL") 445 CASE (4) 446 call wrtout(std_out,"Using interpolation technique mode 1 with full computation of hamiltonian","COLL") 447 CASE DEFAULT 448 write(msg,'(a,i0)')" Wrong interpolation mode for bs_interp_mode: ",dtset%bs_interp_mode 449 MSG_ERROR(msg) 450 END SELECT 451 452 ! Read from dtset 453 if(BSp%use_interp) then 454 BSp%interp_method = dtset%bs_interp_method 455 BSp%rl_nb = dtset%bs_interp_rl_nb 456 BSp%interp_m3_width = dtset%bs_interp_m3_width 457 BSp%interp_kmult(1:3) = dtset%bs_interp_kmult(1:3) 458 BSp%interp_mode = dtset%bs_interp_mode 459 end if 460 461 ! Dimensions and parameters of the calculation. 462 ! TODO one should add npwx as well 463 !BSp%npweps=Dtset%npweps 464 !BSp%npwwfn=Dtset%npwwfn 465 466 ABI_MALLOC(Bsp%lomo_spin, (Bsp%nsppol)) 467 ABI_MALLOC(Bsp%homo_spin, (Bsp%nsppol)) 468 ABI_MALLOC(Bsp%lumo_spin, (Bsp%nsppol)) 469 ABI_MALLOC(Bsp%humo_spin, (Bsp%nsppol)) 470 ABI_MALLOC(Bsp%nbndv_spin, (Bsp%nsppol)) 471 ABI_MALLOC(Bsp%nbndc_spin, (Bsp%nsppol)) 472 473 ! FIXME use bs_loband(nsppol) 474 Bsp%lomo_spin = Dtset%bs_loband 475 write(std_out,*)"bs_loband",Dtset%bs_loband 476 !if (Bsp%nsppol == 2) Bsp%lomo_spin(2) = Dtset%bs_loband 477 478 ! Check lomo correct only for unpolarized semiconductors 479 !if (Dtset%nsppol == 1 .and. Bsp%lomo > Dtset%nelect/2) then 480 ! write(msg,'(a,i0,a,f8.3)') " Bsp%lomo = ",Bsp%lomo," cannot be greater than nelect/2 = ",Dtset%nelect/2 481 ! MSG_ERROR(msg) 482 !end if 483 ! 484 ! ============================================== 485 ! ==== Setup of the q for the optical limit ==== 486 ! ============================================== 487 Bsp%inclvkb = Dtset%inclvkb 488 489 qred2cart = two_pi*Cryst%gprimd 490 qcart2red = qred2cart 491 call matrginv(qcart2red,3,3) 492 493 if (Dtset%gw_nqlwl==0) then 494 BSp%nq = 6 495 ABI_MALLOC(BSp%q,(3,BSp%nq)) 496 BSp%q(:,1) = (/one,zero,zero/) ! (100) 497 BSp%q(:,2) = (/zero,one,zero/) ! (010) 498 BSp%q(:,3) = (/zero,zero,one/) ! (001) 499 BSp%q(:,4) = MATMUL(qcart2red,(/one,zero,zero/)) ! (x) 500 BSp%q(:,5) = MATMUL(qcart2red,(/zero,one,zero/)) ! (y) 501 BSp%q(:,6) = MATMUL(qcart2red,(/zero,zero,one/)) ! (z) 502 else 503 BSp%nq = Dtset%gw_nqlwl 504 ABI_MALLOC(BSp%q,(3,BSp%nq)) 505 BSp%q = Dtset%gw_qlwl 506 end if 507 508 do iq=1,BSp%nq ! normalization 509 qnorm = normv(BSp%q(:,iq),Cryst%gmet,"G") 510 BSp%q(:,iq) = BSp%q(:,iq)/qnorm 511 end do 512 ! 513 ! ====================================================== 514 ! === Define the flags defining the calculation type === 515 ! ====================================================== 516 Bsp%calc_type = Dtset%bs_calctype 517 518 BSp%mbpt_sciss = zero ! Shall we use the scissors operator to open the gap? 519 if (ABS(Dtset%mbpt_sciss)>tol6) BSp%mbpt_sciss = Dtset%mbpt_sciss 520 521 !now test input parameters from input and WFK file and assume some defaults 522 ! 523 ! TODO Add the possibility of using a randomly shifted k-mesh with nsym>1. 524 ! so that densities and potentials are correctly symmetrized but 525 ! the list of the k-point in the IBZ is not expanded. 526 527 if (mband < Dtset%nband(1)) then 528 write(msg,'(2(a,i0),3a,i0)')& 529 & 'WFK file contains only ', mband,' levels instead of ',Dtset%nband(1),' required;',ch10,& 530 & 'The calculation will be done with nbands= ',mband 531 MSG_WARNING(msg) 532 Dtset%nband(:) = mband 533 end if 534 535 BSp%nbnds = Dtset%nband(1) ! TODO Note the change in the meaning of input variables 536 537 if (BSp%nbnds<=Dtset%nelect/2) then 538 write(msg,'(2a,a,i0,a,f8.2)')& 539 & 'BSp%nbnds cannot be smaller than homo ',ch10,& 540 & 'while BSp%nbnds = ',BSp%nbnds,' and Dtset%nelect = ',Dtset%nelect 541 MSG_ERROR(msg) 542 end if 543 544 !TODO add new dim for exchange part and consider the possibility of having npwsigx > npwwfn (see setup_sigma). 545 546 ! === Build enlarged G-sphere for the exchange part === 547 call gsph_extend(Gsph_c,Cryst,Dtset%ecutwfn,Gsph_x) 548 call print_gsphere(Gsph_x,unit=std_out,prtvol=Dtset%prtvol) 549 550 ! NPWVEC as the biggest between npweps and npwwfn. MG RECHECK this part. 551 !BSp%npwwfn = Dtset%npwwfn 552 Bsp%npwwfn = Gsph_x%ng ! FIXME temporary hack 553 BSp%npwvec=MAX(BSp%npwwfn,BSp%npweps) 554 Bsp%ecutwfn = Dtset%ecutwfn 555 556 ! Compute Coulomb term on the largest G-sphere. 557 if (Gsph_x%ng > Gsph_c%ng ) then 558 call vcoul_init(Vcp,Gsph_x,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_x%ng,& 559 & nqlwl,qlwl,ngfftf,comm) 560 else 561 call vcoul_init(Vcp,Gsph_c,Cryst,Qmesh,Kmesh,Dtset%rcut,Dtset%icutcoul,Dtset%vcutgeo,Dtset%ecutsigx,Gsph_c%ng,& 562 & nqlwl,qlwl,ngfftf,comm) 563 end if 564 565 ABI_FREE(qlwl) 566 567 bantot=SUM(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)) 568 ABI_MALLOC(doccde,(bantot)) 569 ABI_MALLOC(eigen,(bantot)) 570 ABI_MALLOC(occfact,(bantot)) 571 doccde=zero; eigen=zero; occfact=zero 572 573 ! Get occupation from input if occopt == 2 574 occ_from_dtset = (Dtset%occopt == 2) 575 576 jj=0; ibtot=0 577 do isppol=1,Dtset%nsppol 578 do ik_ibz=1,Dtset%nkpt 579 do ib=1,Hdr_wfk%nband(ik_ibz+(isppol-1)*Dtset%nkpt) 580 ibtot=ibtot+1 581 if (ib<=BSP%nbnds) then 582 jj=jj+1 583 eigen (jj)=energies_p(ib,ik_ibz,isppol) 584 if (occ_from_dtset) then 585 occfact(jj)=Dtset%occ_orig(ibtot) 586 else 587 occfact(jj)=Hdr_wfk%occ(ibtot) 588 end if 589 end if 590 end do 591 end do 592 end do 593 594 ABI_FREE(energies_p) 595 ! 596 ! * Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of 597 ! symmetry operations in the old GW code (symmorphy and inversion) 598 ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6)) 599 ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt') 600 601 ABI_MALLOC(npwarr,(Dtset%nkpt)) 602 npwarr=BSP%npwwfn 603 604 call ebands_init(bantot,KS_BSt,Dtset%nelect,doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,& 605 & Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,& 606 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 607 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 608 609 ABI_FREE(doccde) 610 ABI_FREE(eigen) 611 ABI_FREE(npwarr) 612 613 !TODO Occupancies are zero if NSCF. One should calculate the occupancies from the energies when 614 ! the occupation scheme for semiconductors is used. 615 call ebands_update_occ(KS_BSt,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 616 617 call ebands_print(KS_BSt,"Band structure read from the WFK file",unit=std_out,prtvol=Dtset%prtvol) 618 619 call ebands_report_gap(KS_BSt,header=" KS band structure",unit=std_out,mode_paral="COLL") 620 621 ABI_MALLOC(val_indeces,(KS_BSt%nkpt,KS_BSt%nsppol)) 622 val_indeces = get_valence_idx(KS_BSt) 623 624 do spin=1,KS_BSt%nsppol 625 val_idx(spin) = val_indeces(1,spin) 626 write(msg,'(a,i2,a,i0)')" For spin : ",spin," val_idx ",val_idx(spin) 627 call wrtout(std_out,msg,"COLL") 628 if ( ANY(val_indeces(1,spin) /= val_indeces(:,spin)) ) then 629 MSG_ERROR("BSE code does not support metals") 630 end if 631 end do 632 633 ABI_FREE(val_indeces) 634 ! 635 ! === Create the BSE header === 636 call hdr_init(KS_BSt,codvsn,Dtset,Hdr_bse,Pawtab,pertcase0,Psps,wvl) 637 638 ! === Get Pawrhoij from the header of the WFK file === 639 ABI_DT_MALLOC(Pawrhoij,(Cryst%natom*Dtset%usepaw)) 640 if (Dtset%usepaw==1) then 641 call pawrhoij_alloc(Pawrhoij,1,Dtset%nspden,Dtset%nspinor,Dtset%nsppol,Cryst%typat,pawtab=Pawtab) 642 call pawrhoij_copy(Hdr_wfk%Pawrhoij,Pawrhoij) 643 end if 644 645 call hdr_update(hdr_bse,bantot,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 646 647 ABI_FREE(occfact) 648 649 if (Dtset%usepaw==1) call pawrhoij_free(Pawrhoij) 650 ABI_DT_FREE(Pawrhoij) 651 652 ! === Find optimal value for G-sphere enlargment due to oscillator matrix elements === 653 654 ! We will split k-points over processors 655 call xmpi_split_work(Kmesh%nbz,comm,my_k1,my_k2,msg,ierr) 656 if (ierr/=0) then 657 MSG_WARNING(msg) 658 end if 659 660 ! If there is no work to do, just skip the computation 661 if (my_k2-my_k1+1 <= 0) then 662 ng0sh_opt(:)=(/zero,zero,zero/) 663 else 664 ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true since bz is buggy 665 ! * -one is used because we loop over all the possibile differences, unlike screening 666 call get_ng0sh(my_k2-my_k1+1,Kmesh%bz(:,my_k1:my_k2),Kmesh%nbz,Kmesh%bz,& 667 & Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt) 668 end if 669 670 call xmpi_max(ng0sh_opt,BSp%mg0,comm,ierr) 671 672 write(msg,'(a,3(i0,1x))') ' optimal value for ng0sh = ',BSp%mg0 673 call wrtout(std_out,msg,"COLL") 674 675 ! === Setup of the FFT mesh for the oscilator strengths === 676 ! * ngfft_osc(7:18)==Dtset%ngfft(7:18) which is initialized before entering screening. 677 ! * Here we redefine ngfft_osc(1:6) according to the following options : 678 ! 679 ! method==0 --> FFT grid read from fft.in (debugging purpose) 680 ! method==1 --> Normal FFT mesh 681 ! method==2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 682 ! method==3 --> Doubled FFT grid, same as the the FFT for the density, 683 ! 684 ! enforce_sym==1 ==> Enforce a FFT mesh compatible with all the symmetry operation and FFT library 685 ! enforce_sym==0 ==> Find the smallest FFT grid compatbile with the library, do not care about symmetries 686 ! 687 ngfft_osc(1:18)=Dtset%ngfft(1:18); method=2 688 if (Dtset%fftgw==00 .or. Dtset%fftgw==01) method=0 689 if (Dtset%fftgw==10 .or. Dtset%fftgw==11) method=1 690 if (Dtset%fftgw==20 .or. Dtset%fftgw==21) method=2 691 if (Dtset%fftgw==30 .or. Dtset%fftgw==31) method=3 692 enforce_sym=MOD(Dtset%fftgw,10) 693 694 call setmesh(gmet,Gsph_x%gvec,ngfft_osc,BSp%npwvec,BSp%npweps,BSp%npwwfn,nfftot_osc,method,BSp%mg0,Cryst,enforce_sym) 695 nfftot_osc=PRODUCT(ngfft_osc(1:3)) 696 697 call print_ngfft(ngfft_osc,"FFT mesh for oscillator matrix elements",std_out,"COLL",prtvol=Dtset%prtvol) 698 ! 699 ! BSp%homo gives the 700 !BSp%homo = val_idx(1) 701 ! highest occupied band for each spin 702 BSp%homo_spin = val_idx 703 704 ! TODO generalize the code to account for this unlikely case. 705 !if (Dtset%nsppol==2) then 706 ! ABI_CHECK(BSp%homo == val_idx(2),"Different valence indeces for spin up and down") 707 !end if 708 709 !BSp%lumo = BSp%homo + 1 710 !BSp%humo = BSp%nbnds 711 !BSp%nbndv = BSp%homo - BSp%lomo + 1 712 !BSp%nbndc = BSp%nbnds - BSp%homo 713 714 BSp%lumo_spin = BSp%homo_spin + 1 715 BSp%humo_spin = BSp%nbnds 716 BSp%nbndv_spin = BSp%homo_spin - BSp%lomo_spin + 1 717 BSp%nbndc_spin = BSp%nbnds - BSp%homo_spin 718 BSp%maxnbndv = MAXVAL(BSp%nbndv_spin(:)) 719 BSp%maxnbndc = MAXVAL(BSp%nbndc_spin(:)) 720 721 BSp%nkbz = Kmesh%nbz 722 723 call ebands_copy(KS_BSt,QP_bst) 724 ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol)) 725 igwene=zero 726 727 call bsp_calctype2str(Bsp,msg) 728 call wrtout(std_out,"Calculation type: "//TRIM(msg)) 729 730 SELECT CASE (Bsp%calc_type) 731 CASE (BSE_HTYPE_RPA_KS) 732 if (ABS(BSp%mbpt_sciss)>tol6) then 733 write(msg,'(a,f8.2,a)')' Applying a scissors operator energy= ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the KS energies." 734 call wrtout(std_out,msg,"COLL") 735 call apply_scissor(QP_BSt,BSp%mbpt_sciss) 736 else 737 write(msg,'(a,f8.2,a)')' Using KS energies since mbpt_sciss= ',BSp%mbpt_sciss*Ha_eV," [eV]." 738 call wrtout(std_out,msg,"COLL") 739 end if 740 741 CASE (BSE_HTYPE_RPA_QPENE) ! Read _GW files with the corrections TODO here I should introduce variable getgw 742 gw_fname=TRIM(Dtfil%filnam_ds(4))//'_GW' 743 gw_fname="__in.gw__" 744 if (.not.file_exists(gw_fname)) then 745 msg = " File "//TRIM(gw_fname)//" not found. Aborting now" 746 MSG_ERROR(msg) 747 end if 748 749 call rdgw(QP_Bst,gw_fname,igwene,extrapolate=.FALSE.) ! here gwenergy is real 750 751 do isppol=1,Dtset%nsppol 752 write(std_out,*) ' k GW energies [eV]' 753 do ik_ibz=1,Kmesh%nibz 754 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(QP_bst%eig(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds) 755 end do 756 write(std_out,*) ' k Im GW energies [eV]' 757 do ik_ibz=1,Kmesh%nibz 758 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(igwene(ib,ik_ibz,isppol)*Ha_eV,ib=1,BSp%nbnds) 759 end do 760 end do 761 ! 762 ! If required apply the scissors operator on top of the QP bands structure (!) 763 if (ABS(BSp%mbpt_sciss)>tol6) then 764 write(msg,'(a,f8.2,a)')' Applying a scissors operator ',BSp%mbpt_sciss*Ha_eV," [eV] on top of the QP energies!" 765 MSG_COMMENT(msg) 766 call apply_scissor(QP_BSt,BSp%mbpt_sciss) 767 end if 768 769 CASE (BSE_HTYPE_RPA_QP) 770 MSG_ERROR("Not implemented error!") 771 772 CASE DEFAULT 773 write(msg,'(a,i0)')"Unknown value for Bsp%calc_type= ",Bsp%calc_type 774 MSG_ERROR(msg) 775 END SELECT 776 777 call ebands_report_gap(QP_BSt,header=" QP band structure",unit=std_out,mode_paral="COLL") 778 779 ! Transitions are ALWAYS ordered in c-v-k mode with k being the slowest index. 780 ! FIXME: linewidths not coded. 781 ABI_MALLOC(gw_energy,(BSp%nbnds,Kmesh%nibz,Dtset%nsppol)) 782 gw_energy = QP_BSt%eig 783 784 BSp%have_complex_ene = ANY(igwene > tol16) 785 786 ! Compute the number of resonant transitions, nreh, for the two spin channels and initialize BSp%Trans. 787 ABI_MALLOC(Bsp%nreh,(Bsp%nsppol)) 788 789 ! Possible cutoff on the transitions. 790 BSp%ircut = Dtset%bs_eh_cutoff(1) 791 BSp%uvcut = Dtset%bs_eh_cutoff(2) 792 793 call init_transitions(BSp%Trans,BSp%lomo_spin,BSp%humo_spin,BSp%ircut,Bsp%uvcut,BSp%nkbz,Bsp%nbnds,Bsp%nkibz,& 794 & BSp%nsppol,Dtset%nspinor,gw_energy,QP_BSt%occ,Kmesh%tab,minmax_tene,Bsp%nreh) 795 796 ! Setup of the frequency mesh for the absorption spectrum. 797 ! If not specified, use the min-max resonant transition energy and make it 10% smaller|larger. 798 799 !if (ABS(Dtset%bs_freq_mesh(1)) < tol6) then 800 ! Dtset%bs_freq_mesh(1) = MAX(minmax_tene(1) - minmax_tene(1) * 0.1, zero) 801 !end if 802 803 if (ABS(Dtset%bs_freq_mesh(2)) < tol6) then 804 Dtset%bs_freq_mesh(2) = minmax_tene(2) + minmax_tene(2) * 0.1 805 end if 806 807 Bsp%omegai = Dtset%bs_freq_mesh(1) 808 Bsp%omegae = Dtset%bs_freq_mesh(2) 809 Bsp%domega = Dtset%bs_freq_mesh(3) 810 BSp%broad = Dtset%zcut 811 812 ! The frequency mesh (including the complex imaginary shift) 813 BSp%nomega = (BSp%omegae - BSp%omegai)/BSp%domega + 1 814 ABI_MALLOC(BSp%omega,(BSp%nomega)) 815 do io=1,BSp%nomega 816 BSp%omega(io) = (BSp%omegai + (io-1)*BSp%domega) + j_dpc*BSp%broad 817 end do 818 819 ABI_FREE(gw_energy) 820 ABI_FREE(igwene) 821 822 do spin=1,Bsp%nsppol 823 write(msg,'(a,i2,a,i0)')" For spin: ",spin,' the number of resonant e-h transitions is: ',BSp%nreh(spin) 824 call wrtout(std_out,msg,"COLL") 825 end do 826 827 if (ANY(Bsp%nreh/=Bsp%nreh(1))) then 828 write(msg,'(a,2(i0,1x))')" BSE code with different number of transitions for the two spin channels: ",Bsp%nreh 829 MSG_WARNING(msg) 830 end if 831 ! 832 ! Create transition table vcks2t 833 Bsp%lomo_min = MINVAL(BSp%lomo_spin) 834 Bsp%homo_max = MAXVAL(BSp%homo_spin) 835 Bsp%lumo_min = MINVAL(BSp%lumo_spin) 836 Bsp%humo_max = MAXVAL(BSp%humo_spin) 837 838 ABI_MALLOC(Bsp%vcks2t,(BSp%lomo_min:BSp%homo_max,BSp%lumo_min:BSp%humo_max,BSp%nkbz,Dtset%nsppol)) 839 Bsp%vcks2t = 0 840 841 do spin=1,BSp%nsppol 842 do it=1,BSp%nreh(spin) 843 BSp%vcks2t(BSp%Trans(it,spin)%v,BSp%Trans(it,spin)%c,BSp%Trans(it,spin)%k,spin) = it 844 end do 845 end do 846 847 hexc_size = SUM(Bsp%nreh); if (Bsp%use_coupling>0) hexc_size=2*hexc_size 848 if (Bsp%nstates<=0) then 849 Bsp%nstates=hexc_size 850 else 851 if (Bsp%nstates>hexc_size) then 852 Bsp%nstates=hexc_size 853 write(msg,'(2(a,i0),2a)')& 854 & "Since the total size of excitonic Hamiltonian ",hexc_size," is smaller than Bsp%nstates ",Bsp%nstates,ch10,& 855 & "the number of excitonic states nstates has been modified" 856 MSG_WARNING(msg) 857 end if 858 end if 859 860 msg=' Fundamental parameters for the solution of the Bethe-Salpeter equation:' 861 call print_bs_parameters(BSp,unit=std_out,header=msg,mode_paral="COLL",prtvol=Dtset%prtvol) 862 call print_bs_parameters(BSp,unit=ab_out, header=msg,mode_paral="COLL") 863 864 if (ANY (Cryst%symrec(:,:,1) /= RESHAPE ( (/1,0,0,0,1,0,0,0,1/),(/3,3/) )) .or. & 865 & ANY( ABS(Cryst%tnons(:,1)) > tol6) ) then 866 write(msg,'(3a,9i2,2a,3f6.3,2a)')& 867 & "The first symmetry operation should be the Identity with zero tnons while ",ch10,& 868 & "symrec(:,:,1) = ",Cryst%symrec(:,:,1),ch10,& 869 & "tnons(:,1) = ",Cryst%tnons(:,1),ch10,& 870 & "This is not allowed, sym_rhotwgq0 should be changed." 871 MSG_ERROR(msg) 872 end if 873 ! 874 ! Prefix for generic output files. 875 BS_files%out_basename = TRIM(Dtfil%filnam_ds(4)) 876 ! 877 ! Search for files to restart from. 878 if (Dtset%gethaydock/=0 .or. Dtset%irdhaydock/=0) then 879 BS_files%in_haydock_basename = TRIM(Dtfil%fnameabi_haydock) 880 end if 881 882 test_file = Dtfil%fnameabi_bsham_reso 883 if (file_exists(test_file)) then 884 BS_files%in_hreso = test_file 885 else 886 BS_files%out_hreso = TRIM(Dtfil%filnam_ds(4))//'_BSR' 887 end if 888 889 test_file = Dtfil%fnameabi_bsham_coup 890 if (file_exists(test_file) ) then 891 BS_files%in_hcoup = test_file 892 else 893 BS_files%out_hcoup = TRIM(Dtfil%filnam_ds(4))//'_BSC' 894 end if 895 ! 896 ! in_eig is the name of the input file with eigenvalues and eigenvectors 897 ! constructed from getbseig or irdbseig. out_eig is the name of the output file 898 ! produced by this dataset. in_eig_exists checks for the presence of the input file. 899 ! 900 if (file_exists(Dtfil%fnameabi_bseig)) then 901 BS_files%in_eig = Dtfil%fnameabi_bseig 902 else 903 BS_files%out_eig = TRIM(BS_files%out_basename)//"_BSEIG" 904 end if 905 906 call print_bs_files(BS_files,unit=std_out) 907 908 909 ! 910 ! ========================================================== 911 ! ==== Temperature dependence of the spectrum ============== 912 ! ========================================================== 913 BSp%do_ep_renorm = .FALSE. 914 BSp%do_lifetime = .FALSE. ! Not yet implemented 915 916 ep_nc_fname = 'test_EP.nc' 917 if(file_exists(ep_nc_fname)) then 918 BSp%do_ep_renorm = .TRUE. 919 920 if(my_rank == master) then 921 call eprenorms_from_epnc(Epren,ep_nc_fname) 922 end if 923 call eprenorms_bcast(Epren,master,comm) 924 end if 925 926 927 ! 928 ! ========================================================== 929 ! ==== Final check on the parameters of the calculation ==== 930 ! ========================================================== 931 if ( Bsp%use_coupling>0 .and. ALL(Bsp%algorithm /= [BSE_ALGO_DDIAGO, BSE_ALGO_HAYDOCK]) ) then 932 msg = "Resonant+Coupling is only available with the direct diagonalization or the haydock method." 933 MSG_ERROR(msg) 934 end if 935 936 ! autoparal section 937 if (dtset%max_ncpus /=0 .and. dtset%autoparal /=0 ) then 938 ount = ab_out 939 ! TODO: 940 ! nsppol and calculation with coupling! 941 942 ! Temporary table needed to estimate memory 943 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 944 if (Dtset%usepaw==1) then 945 do iat=1,Cryst%natom 946 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 947 end do 948 end if 949 950 tot_nreh = SUM(BSp%nreh) 951 work_size = tot_nreh * (tot_nreh + 1) / 2 952 953 write(ount,'(a)')"--- !Autoparal" 954 write(ount,"(a)")'#Autoparal section for Bethe-Salpeter runs.' 955 956 write(ount,"(a)") "info:" 957 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 958 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 959 write(ount,"(a,i0)")" nkibz: ",Bsp%nkibz 960 write(ount,"(a,i0)")" nkbz: ",Bsp%nkbz 961 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 962 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 963 write(ount,"(a,i0)")" lomo_min: ",Bsp%lomo_min 964 write(ount,"(a,i0)")" humo_max: ",Bsp%humo_max 965 write(ount,"(a,i0)")" tot_nreh: ",tot_nreh 966 !write(ount,"(a,i0)")" nbnds: ",Ep%nbnds 967 968 ! Wavefunctions are not distributed. We read all the bands 969 ! from 1 up to Bsp%nbnds because we have to recompute rhor 970 ! but then we deallocate all the states that are not used for the construction of the e-h 971 ! before allocating the EXC hamiltonian. Hence we can safely use (humo - lomo + 1) instead of Bsp%nbnds. 972 !my_nbks = (Bsp%humo - Bsp%lomo +1) * Bsp%nkibz * Dtset%nsppol 973 974 ! This one overestimates the memory but it seems to be safer. 975 my_nbks = Bsp%nbnds * Dtset%nkpt * Dtset%nsppol 976 977 ! Memory needed for Fourier components ug. 978 ug_mem = two*gwpc*Dtset%nspinor*Bsp%npwwfn*my_nbks*b2Mb 979 980 ! Memory needed for real space ur. 981 ur_mem = zero 982 if (MODULO(Dtset%gwmem,10)==1) then 983 ur_mem = two*gwpc*Dtset%nspinor*nfftot_osc*my_nbks*b2Mb 984 end if 985 986 ! Memory needed for PAW projections Cprj 987 cprj_mem = zero 988 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 989 990 wfsmem_mb = ug_mem + ur_mem + cprj_mem 991 992 ! Non-scalable memory in Mb i.e. memory that is not distributed with MPI: wavefunctions + W 993 nonscal_mem = (wfsmem_mb + two*gwpc*BSp%npweps**2*b2Mb) * 1.1_dp 994 995 ! List of configurations. 996 write(ount,"(a)")"configurations:" 997 do il=1,dtset%max_ncpus 998 if (il > work_size) cycle 999 neh_per_proc = work_size / il 1000 neh_per_proc = neh_per_proc + MOD(work_size, il) 1001 eff = (one * work_size) / (il * neh_per_proc) 1002 1003 ! EXC matrix is distributed. 1004 mempercpu_mb = nonscal_mem + two * dpc * neh_per_proc * b2Mb 1005 1006 write(ount,"(a,i0)")" - tot_ncpus: ",il 1007 write(ount,"(a,i0)")" mpi_ncpus: ",il 1008 !write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 1009 write(ount,"(a,f12.9)")" efficiency: ",eff 1010 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 1011 end do 1012 1013 write(ount,'(a)')"..." 1014 1015 ABI_FREE(nlmn_atm) 1016 MSG_ERROR_NODUMP("aborting now") 1017 end if 1018 1019 DBG_EXIT("COLL") 1020 1021 end subroutine setup_bse