TABLE OF CONTENTS
- ABINIT/m_sigma_driver
- ABINIT/paw_qpscgw
- ABINIT/setup_sigma
- ABINIT/setup_vcp
- ABINIT/sigma_bksmask
- ABINIT/sigma_tables
- m_sigma_driver/sigma
ABINIT/m_sigma_driver [ Modules ]
NAME
m_sigma_driver
FUNCTION
Calculate the matrix elements of the self-energy operator.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (MG, GMR, VO, LR, RWG, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_sigma_driver 23 24 use defs_basis 25 use m_gwdefs 26 use defs_wvltypes 27 use m_dtset 28 use m_xmpi 29 use m_xomp 30 use m_errors 31 use m_abicore 32 use m_ab7_mixing 33 use m_nctk 34 use m_kxc 35 use m_distribfft 36 use netcdf 37 use m_hdr 38 use libxc_functionals 39 use m_dtfil 40 use m_crystal 41 use m_cgtools 42 43 use defs_datatypes, only : pseudopotential_type, ebands_t 44 use defs_abitypes, only : MPI_type 45 use m_time, only : timab 46 use m_numeric_tools, only : imax_loc 47 use m_fstrings, only : strcat, sjoin, itoa, basename, ktoa, ltoa 48 use m_hide_blas, only : xdotc 49 use m_io_tools, only : open_file, file_exists, iomode_from_fname 50 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 51 use m_geometry, only : normv, mkrdim, metric 52 use m_fftcore, only : print_ngfft 53 use m_fft_mesh, only : get_gftt, setmesh 54 use m_fft, only : fourdp 55 use m_ioarr, only : fftdatar_write, read_rhor 56 use m_ebands, only : ebands_update_occ, ebands_copy, ebands_report_gap, ebands_get_valence_idx, ebands_get_bandenergy,& 57 ebands_free, ebands_init, ebands_ncwrite, ebands_interpolate_kpath, get_eneocc_vect, & 58 ebands_enclose_degbands, ebands_get_gaps, gaps_t 59 use m_energies, only : energies_type, energies_init 60 use m_bz_mesh, only : kmesh_t, littlegroup_t, littlegroup_free, isamek, get_ng0sh, find_qmesh 61 use m_gsphere, only : gsphere_t, merge_and_sort_kg, setshells 62 use m_kg, only : getph, getcut 63 use m_xcdata, only : get_xclevel 64 use m_wfd, only : wfd_init, wfdgw_t, wfdgw_copy, test_charge, wave_t 65 use m_vcoul, only : vcoul_t 66 use m_qparticles, only : wrqps, rdqps, rdgw, show_QP, updt_m_ks_to_qp 67 use m_screening, only : mkdump_er, em1results_free, epsilonm1_results, init_er_from_file 68 use m_ppmodel, only : ppm_init, ppm_free, setup_ppmodel, getem1_from_PPm, ppmodel_t 69 use m_sigma, only : sigma_init, sigma_free, sigma_ncwrite, sigma_t, sigma_get_exene, & 70 mels_get_haene, mels_get_kiene, sigma_get_excene, write_sigma_header, write_sigma_results 71 use m_dyson_solver, only : solve_dyson 72 use m_esymm, only : esymm_t, esymm_free, esymm_failed 73 use m_melemts, only : melflags_t, melements_t 74 use m_pawang, only : pawang_type 75 use m_pawrad, only : pawrad_type 76 use m_pawtab, only : pawtab_type, pawtab_print, pawtab_get_lsize 77 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 78 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print 79 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print 80 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_copy, pawrhoij_free, & 81 pawrhoij_inquire_dim, pawrhoij_symrhoij, pawrhoij_unpack 82 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, paw_overlap 83 use m_pawdij, only : pawdij, symdij_all 84 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 85 use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 86 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free 87 use m_paw_slater, only : paw_mkdijexc_core, paw_dijhf 88 use m_paw_dmft, only : paw_dmft_type 89 use m_paw_sphharm, only : setsym_ylm 90 use m_paw_mkrho, only : denfgr 91 use m_paw_nhat, only : nhatgrid, pawmknhat 92 use m_paw_tools, only : chkpawovlp, pawprt 93 use m_paw_denpot, only : pawdenpot 94 use m_paw_init, only : pawinit, paw_gencond 95 use m_classify_bands,only : classify_bands 96 use m_wfk, only : wfk_read_eigenvalues 97 use m_io_kss, only : make_gvec_kss 98 use m_vhxc_me, only : calc_vhxc_me 99 use m_cohsex, only : cohsex_me 100 use m_sigx, only : calc_sigx_me 101 use m_sigc, only : calc_sigc_me 102 use m_setvtr, only : setvtr 103 use m_mkrho, only : prtrhomxmn 104 use m_pspini, only : pspini 105 use m_calc_ucrpa, only : calc_ucrpa 106 use m_prep_calc_ucrpa,only : prep_calc_ucrpa 107 use m_paw_correlations,only : pawpuxinit 108 use m_spacepar, only : hartre 109 use m_gwrdm, only : calc_rdmx,calc_rdmc,natoccs,update_hdr_bst,print_tot_occ,get_chkprdm,& 110 print_chkprdm,change_matrix,print_total_energy,print_band_energies,quadrature_sigma_cw 111 use m_plowannier, only : operwan_realspace_type,plowannier_type,init_plowannier,get_plowannier,& 112 fullbz_plowannier,init_operwan_realspace,reduce_operwan_realspace,& 113 destroy_operwan_realspace,destroy_plowannier,zero_operwan_realspace 114 115 implicit none 116 117 private
ABINIT/paw_qpscgw [ Functions ]
NAME
paw_qpscgw
FUNCTION
This routine is called during QP self-consistent GW calculations. It calculates the new QP on-site quantities using the QP amplitudes read from the QPS file.
INPUTS
Wfd<wfdgw_t>=Datatype gathering data on QP amplitudes. nscf=Number of QPSCGW iterations done so far (read from the QPS file). nfftf=Number of points in the fine FFT grid. ngfft(18)=information on the fine FFT grid used for densities and potentials. Dtset<dataset_type>=All input variables for this dataset. Cryst<crystal_t>=Info on unit cell and symmetries. Kmesh<kmesh_t>=Structure describing the k-point sampling. Psps<Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the KSS file Pawang<pawang_type>=PAW angular mesh and related data. Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data Pawfgrtab(natom)<Pawfgrtab_type>= For PAW, various arrays giving data related to fine grid for a given atom. prev_Pawrhoij(Cryst%natom))<Pawrhoij_type>=Previous QP rhoij used for mixing if nscf>0 and rhoqpmix/=one. MPI_enreg=Information about MPI parallelization. qp_ebands<ebands_t>=QP band structure. QP_energies<Energies_type>=Simple datastructure to gather all part of total energy. nhatgrdim= 0 if pawgrnhat array is not used ; 1 otherwise
OUTPUT
QP_pawrhoij(Cryst%natom))<Pawrhoij_type>=on-site densities calculated from the QP amplitudes. qp_nhat(nfftf,Dtset%nspden)=Compensation charge density calculated from the QP amplitudes. qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim)=Derivatives of the QP nhat on fine rectangular grid (and derivatives). qp_compch_sph=QP compensation charge integral inside spheres computed over spherical meshes. qp_compch_fft=QP compensation charge inside spheres computed over fine fft grid. QP_paw_ij(Cryst%natom)<Paw_ij_type>=Non-local D_ij strengths of the QP Hamiltonian. QP_paw_an(Cryst%natom)<Paw_an_type>=Various arrays related to the Hamiltonian given on ANgular mesh or ANgular moments.
SOURCE
4161 subroutine paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,qp_ebands, & 4162 Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij, & 4163 QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim, & 4164 qp_nhatgr,qp_compch_sph,qp_compch_fft) 4165 4166 !Arguments ------------------------------------ 4167 !scalars 4168 integer,intent(in) :: nfftf,nscf,nhatgrdim 4169 real(dp),intent(out) :: qp_compch_fft,qp_compch_sph 4170 type(kmesh_t),intent(in) :: Kmesh 4171 type(crystal_t),intent(in) :: Cryst 4172 type(Dataset_type),intent(in) :: Dtset 4173 type(Pseudopotential_type),intent(in) :: Psps 4174 type(Pawang_type),intent(in) :: Pawang 4175 type(ebands_t),intent(in) :: qp_ebands 4176 type(Energies_type),intent(inout) :: QP_energies 4177 type(wfdgw_t),intent(inout) :: Wfd 4178 !arrays 4179 integer,intent(in) :: ngfftf(18) 4180 real(dp),intent(out) :: qp_nhat(nfftf,Dtset%nspden) 4181 real(dp),intent(out) :: qp_nhatgr(nfftf,Dtset%nspden,3*nhatgrdim) 4182 type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom) 4183 type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw) 4184 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 4185 type(Pawrhoij_type),intent(inout) :: prev_Pawrhoij(Cryst%natom) 4186 type(Pawrhoij_type),intent(out) :: QP_pawrhoij(Cryst%natom) 4187 type(Paw_ij_type),intent(out) :: QP_paw_ij(Cryst%natom) 4188 type(Paw_an_type),intent(inout) :: QP_paw_an(Cryst%natom) 4189 4190 !Local variables ------------------------------ 4191 !scalars 4192 integer,parameter :: ipert0 = 0, idir0 = 0, optrhoij1 = 1 4193 integer :: choice,cplex,cplex_rhoij,has_dijU,has_dijso,iat,ider 4194 integer :: izero,nkxc1,nspden_rhoij,nzlmopt, option,usexcnhat 4195 character(len=500) :: msg 4196 !arrays 4197 real(dp) :: k0(3) 4198 4199 !************************************************************************ 4200 4201 ABI_UNUSED(Kmesh%nibz) 4202 ! 4203 ! * 0 if Vloc in atomic data is Vbare (Blochl s formulation) 4204 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse s formulation) 4205 usexcnhat=MAXVAL(Pawtab(:)%usexcnhat) 4206 ! 4207 ! Calculate new rhoij_qp from updated Cprj_ibz, note use_rhoij_=1. 4208 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 4209 nspden=Dtset%nspden,spnorb=Dtset%pawspnorb,cpxocc=Dtset%pawcpxocc) 4210 call pawrhoij_alloc(QP_pawrhoij,cplex_rhoij,nspden_rhoij,Dtset%nspinor,Dtset%nsppol,Cryst%typat,& 4211 pawtab=Pawtab,use_rhoij_=1,use_rhoijres=1) 4212 4213 ! FIXME kptop should be passed via Kmesh, in GW time reversal is always assumed. 4214 call wfd%pawrhoij(Cryst,qp_ebands,Dtset%kptopt,QP_pawrhoij,Dtset%pawprtvol) 4215 ! 4216 ! Symmetrize QP $\rho_{ij}$. 4217 choice=1 4218 call pawrhoij_symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,& 4219 Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij1,Pawang,Dtset%pawprtvol,& 4220 Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 4221 4222 ! ====================== 4223 ! ==== Make QP nhat ==== 4224 ! ====================== 4225 cplex=1; ider=2*nhatgrdim; izero=0; k0(:)=zero 4226 4227 call pawmknhat(qp_compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,& 4228 Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 4229 Pawfgrtab,qp_nhatgr,qp_nhat,QP_Pawrhoij,QP_Pawrhoij,Pawtab,k0,Cryst%rprimd,& 4230 Cryst%ucvol,Dtset%usewvl,Cryst%xred) 4231 4232 ! Allocate quantities related to the PAW spheres for the QP Hamiltonian. 4233 ! TODO call paw_ij_init in scfcv and respfn, fix small issues 4234 has_dijso=Dtset%pawspnorb; has_dijU=merge(0,1,Dtset%usepawu==0) 4235 4236 call paw_ij_nullify(QP_paw_ij) 4237 call paw_ij_init(QP_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,& 4238 Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 4239 has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,& 4240 has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1) 4241 4242 call paw_an_nullify(QP_paw_an); nkxc1=0 ! No kernel 4243 call paw_an_init(QP_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,& 4244 cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1) 4245 4246 ! ===================================================== 4247 ! ==== Optional mixing of the PAW onsite densities ==== 4248 ! ===================================================== 4249 if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then 4250 write(msg,'(a,f6.3)')' sigma: mixing on-site QP rho_ij densities using rhoqpmix= ',Dtset%rhoqpmix 4251 call wrtout(std_out, msg) 4252 ! qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor) 4253 4254 call pawrhoij_unpack(QP_Pawrhoij) ! Unpack new QP %rhoijp 4255 call pawrhoij_unpack(prev_Pawrhoij) ! Unpack previous QP %rhoijp 4256 4257 do iat=1,Cryst%natom 4258 QP_pawrhoij(iat)%rhoij_ = prev_Pawrhoij(iat)%rhoij_ & 4259 + Dtset%rhoqpmix * (QP_pawrhoij(iat)%rhoij_ - prev_pawrhoij(iat)%rhoij_) 4260 4261 prev_pawrhoij(iat)%use_rhoij_=0 4262 ABI_FREE(prev_pawrhoij(iat)%rhoij_) 4263 end do 4264 4265 ! Re-Symmetrize mixed QP $\rho_{ij}$. 4266 choice=1 4267 call pawrhoij_symrhoij(QP_pawrhoij,QP_pawrhoij,choice,Cryst%gprimd,Cryst%indsym,ipert0,& 4268 Cryst%natom,Cryst%nsym,Cryst%ntypat,optrhoij1,Pawang,Dtset%pawprtvol,& 4269 Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,Cryst%typat) 4270 end if 4271 4272 do iat=1,Cryst%natom 4273 QP_pawrhoij(iat)%use_rhoij_=0 4274 ABI_FREE(QP_pawrhoij(iat)%rhoij_) 4275 end do 4276 ! 4277 ! ================================================================================= 4278 ! ==== Evaluate on-site energies, potentials, densities using (mixed) QP rhoij ==== 4279 ! ================================================================================= 4280 ! * Initialize also "lmselect" (index of non-zero LM-moments of densities). 4281 4282 nzlmopt=-1; option=0; qp_compch_sph=greatest_real 4283 4284 call pawdenpot(qp_compch_sph,QP_energies%e_paw,QP_energies%e_pawdc,& 4285 ipert0,Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 4286 Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,QP_paw_an,QP_paw_an,& 4287 QP_paw_ij,Pawang,Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,& 4288 Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,& 4289 Cryst%ucvol,Psps%znuclpsp) 4290 4291 end subroutine paw_qpscgw
ABINIT/setup_sigma [ Functions ]
NAME
setup_sigma
FUNCTION
Initialize the data type containing parameters for a sigma calculation.
INPUTS
acell(3)=length scales of primitive translations (bohr) wfk_fname=Name of the WFK file. Dtset<type(dataset_type)>=all input variables for this dataset Dtfil<type(datafiles_type)>=variables related to files rprim(3,3)=dimensionless real space primitive translations Psps <Pseudopotential_type)>=Info on pseudopotential, only for consistency check of the WFK file
OUTPUT
Sigp<sigparams_t>=Parameters governing the self-energy calculation. Kmesh <kmesh_t>=Structure describing the k-point sampling. Qmesh <kmesh_t>=Structure describing the q-point sampling. Cryst<crystal_t>=Info on unit cell and symmetries. Gsph_Max<gsphere_t>=Info on the G-sphere Gsph_c<gsphere_t>=Info on the G-sphere for W and Sigma_c Gsph_x<gsphere_t>=Info on the G-sphere for and Sigma_x Hdr_wfk<hdr_type>=The header of the WFK file Hdr_out<hdr_type>=The header to be used for the results of sigma calculations. Vcp<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique. Er<Epsilonm1_results>=Datatype storing data used to construct the screening (partially Initialized in OUTPUT) ks_ebands<ebands_t>=The KS energies and occupation factors. gwc_ngfft(18), gwx_ngfft(18)= FFT meshes for the oscillator strengths used for the correlated and the exchange part of the self-energy, respectively. comm=MPI communicator.
SOURCE
2990 subroutine setup_sigma(codvsn,wfk_fname,acell,rprim,Dtset,Dtfil,Psps,Pawtab,& 2991 & gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_out,Cryst,Kmesh,Qmesh,ks_ebands,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm) 2992 2993 !Arguments ------------------------------------ 2994 !scalars 2995 integer,intent(in) :: comm 2996 character(len=8),intent(in) :: codvsn 2997 character(len=*),intent(in) :: wfk_fname 2998 type(Datafiles_type),intent(in) :: Dtfil 2999 type(Dataset_type),intent(inout) :: Dtset 3000 type(Pseudopotential_type),intent(in) :: Psps 3001 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Dtset%usepaw) 3002 type(sigparams_t),intent(out) :: Sigp 3003 type(Epsilonm1_results),intent(out) :: Er 3004 type(ebands_t),intent(out) :: ks_ebands 3005 type(kmesh_t),intent(out) :: Kmesh,Qmesh 3006 type(crystal_t),intent(out) :: Cryst 3007 type(gsphere_t),intent(out) :: Gsph_Max,Gsph_x,Gsph_c 3008 type(Hdr_type),intent(out) :: Hdr_wfk,Hdr_out 3009 type(vcoul_t),intent(out) :: Vcp 3010 !arrays 3011 integer,intent(out) :: gwc_ngfft(18),gwx_ngfft(18) 3012 real(dp),intent(in) :: acell(3),rprim(3,3) 3013 3014 !Local variables------------------------------- 3015 !scalars 3016 integer,parameter :: pertcase0=0,master=0 3017 integer :: bantot,enforce_sym,gwcalctyp,ib,ibtot,icutcoul_eff,ii,ikcalc,ikibz,io,isppol,itypat,jj,method 3018 integer :: mod10,mqmem,mband,ng_kss,nsheps,ikcalc2bz,ierr,gap_err,ng, nsppol 3019 integer :: gwc_nfftot,gwx_nfftot,nqlwl,test_npwkss,my_rank,nprocs,ik,nk_found,ifo,timrev,usefock_ixc 3020 integer :: iqbz,isym,iq_ibz,itim,ic,pinv,ig1,ng_sigx,spin,gw_qprange,ivcoul_init,nvcoul_init,xclevel_ixc 3021 real(dp),parameter :: OMEGASIMIN=0.01d0,tol_enediff=0.001_dp*eV_Ha 3022 real(dp) :: domegas,domegasi,ucvol,rcut 3023 logical,parameter :: linear_imag_mesh=.TRUE. 3024 logical :: ltest,remove_inv,changed,found 3025 character(len=500) :: msg 3026 character(len=fnlen) :: fname,fcore,string 3027 type(wvl_internal_type) :: wvl 3028 type(gaps_t) :: gaps 3029 !arrays 3030 integer :: ng0sh_opt(3),G0(3),q_umklp(3),kpos(6), units(2) 3031 integer,allocatable :: npwarr(:),val_indices(:,:) 3032 integer,pointer :: gvec_kss(:,:),gsphere_sigx_p(:,:) 3033 integer,pointer :: test_gvec_kss(:,:) 3034 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),sq(3),q_bz(3),gamma_point(3,1) 3035 real(dp),pointer :: energies_p(:,:,:) 3036 real(dp),allocatable :: doccde(:),eigen(:),occfact(:),qlwl(:,:) 3037 type(Pawrhoij_type),allocatable :: Pawrhoij(:) 3038 type(vcoul_t) :: Vcp_ks 3039 3040 ! ************************************************************************* 3041 3042 DBG_ENTER('COLL') 3043 units = [std_out, ab_out] 3044 3045 ! Check for calculations that are not implemented 3046 nsppol = dtset%nsppol 3047 ltest = ALL(Dtset%nband(1:Dtset%nkpt*nsppol) == Dtset%nband(1)) 3048 ABI_CHECK(ltest,'Dtset%nband(:) must be constant') 3049 3050 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 3051 3052 ! Basic parameters 3053 Sigp%ppmodel = Dtset%ppmodel 3054 Sigp%gwcalctyp = Dtset%gwcalctyp 3055 Sigp%nbnds = Dtset%nband(1) 3056 Sigp%symsigma = Dtset%symsigma 3057 Sigp%zcut = Dtset%zcut 3058 Sigp%mbpt_sciss = Dtset%mbpt_sciss 3059 3060 timrev= 2 ! This information is not reported in the header 3061 ! 1 => do not use time-reversal symmetry 3062 ! 2 => take advantage of time-reversal symmetry 3063 if (any(dtset%kptopt == [3, 4])) timrev = 1 3064 3065 ! === For HF, SEX or COHSEX use Hybertsen-Louie PPM (only $\omega=0$) === 3066 ! * Use fake screening for HF. 3067 ! FIXME Why, we should not redefine Sigp%ppmodel 3068 gwcalctyp=Sigp%gwcalctyp 3069 mod10 =MOD(Sigp%gwcalctyp,10) 3070 if (mod10==5.or.mod10==6.or.mod10==7) Sigp%ppmodel=2 3071 if (mod10<5.and.MOD(Sigp%gwcalctyp,1)/=1) then ! * One shot GW (PPM or contour deformation). 3072 if (Dtset%nomegasrd==1) then ! avoid division by zero! 3073 Sigp%nomegasrd =1 3074 Sigp%maxomega4sd=zero 3075 Sigp%deltae =zero 3076 else 3077 Sigp%nomegasrd = Dtset%nomegasrd 3078 Sigp%maxomega4sd = Dtset%omegasrdmax 3079 Sigp%deltae = (2*Sigp%maxomega4sd)/(Sigp%nomegasrd-1) 3080 endif 3081 else 3082 ! For AC no need to evaluate derivative by finite differences. 3083 Sigp%nomegasrd =1 3084 Sigp%maxomega4sd=zero 3085 Sigp%deltae =zero 3086 end if 3087 3088 ! For analytic continuation define the number of imaginary frequencies for Sigma 3089 ! Tests show than more than 12 freqs in the Pade approximant worsen the results! 3090 Sigp%nomegasi=0 3091 3092 if (mod10==1) then 3093 Sigp%nomegasi =Dtset%nomegasi 3094 Sigp%omegasimax=Dtset%omegasimax 3095 Sigp%omegasimin=OMEGASIMIN 3096 write(msg,'(4a,i3,2(2a,f8.3),a)')ch10,& 3097 ' Parameters for analytic continuation : ',ch10,& 3098 ' number of imaginary frequencies for sigma = ',Sigp%nomegasi,ch10,& 3099 ' min frequency for sigma on imag axis [eV] = ',Sigp%omegasimin*Ha_eV,ch10,& 3100 ' max frequency for sigma on imag axis [eV] = ',Sigp%omegasimax*Ha_eV,ch10 3101 call wrtout(std_out, msg) 3102 3103 !TODO this should not be done here but in init_sigma_t 3104 ABI_MALLOC(Sigp%omegasi,(Sigp%nomegasi)) 3105 3106 if (linear_imag_mesh) then 3107 ! Linear mesh along the imaginary axis. 3108 domegasi=Sigp%omegasimax/(Sigp%nomegasi-1) 3109 do io=1,Sigp%nomegasi 3110 Sigp%omegasi(io)=CMPLX(zero,(io-1)*domegasi) 3111 end do 3112 else 3113 ! Logarithmic mesh along the imaginary axis. 3114 ABI_ERROR("AC + log mesh not implemented") 3115 !domegasi=(Sigp%omegasimax/Sigp%omegasimin)**(one/(Sigp%nomegasi-1)) 3116 !Sigp%omegasi(1)=czero; ldi=domegasi 3117 !do io=2,Sigp%nomegasi 3118 ! omega(io)=CMPLX(zero,ldi*Sigp%omegasimin) 3119 ! Sigp%omegasi(io)=ldi*domegasi 3120 !end do 3121 end if 3122 3123 ! MRM: do not print for 1-RDM correction 3124 if(Sigp%gwcalctyp/=21) then 3125 write(msg,'(4a)')ch10,& 3126 ' setup_sigma : calculating Sigma(iw)',& 3127 ' at imaginary frequencies [eV] (Fermi Level set to 0) ',ch10 3128 call wrtout(units, msg) 3129 do io=1,Sigp%nomegasi 3130 write(msg,'(2(f10.3,2x))')Sigp%omegasi(io)*Ha_eV 3131 call wrtout(units, msg) 3132 end do 3133 endif 3134 3135 ltest=(Sigp%omegasimax>0.1d-4.and.Sigp%nomegasi>0) 3136 ABI_CHECK(ltest,'Wrong value of omegasimax or nomegasi') 3137 if (Sigp%gwcalctyp/=1) then ! only one shot GW is allowed for AC. 3138 !ABI_ERROR("SC-GW with analytic continuation is not coded") ! MRM: let's allow it 3139 end if 3140 end if 3141 3142 if (Sigp%symsigma/=0.and.gwcalctyp>=20) then 3143 ABI_WARNING("SC-GW with symmetries is still under development. Use at your own risk!") 3144 ABI_ERROR("SC-GW requires symsigma == 0 in input. New default in Abinit9 is symsigma 1!") 3145 end if 3146 3147 ! Setup parameters for Spectral function. 3148 if (Dtset%gw_customnfreqsp/=0) then 3149 Sigp%nomegasr = Dtset%gw_customnfreqsp 3150 ABI_WARNING('Custom grid for spectral function specified. Assuming experienced user.') 3151 if (Dtset%gw_customnfreqsp/=0) then 3152 Dtset%nfreqsp = Dtset%gw_customnfreqsp 3153 ABI_WARNING('nfreqsp has been set to the same number as gw_customnfreqsp') 3154 end if 3155 else 3156 Sigp%nomegasr =Dtset%nfreqsp 3157 Sigp%minomega_r=Dtset%freqspmin 3158 Sigp%maxomega_r=Dtset%freqspmax 3159 end if 3160 3161 if (Sigp%nomegasr>0) then 3162 if (Dtset%gw_customnfreqsp==0) then 3163 ! Check 3164 if (Sigp%minomega_r >= Sigp%maxomega_r) then 3165 ABI_ERROR('freqspmin must be smaller than freqspmax!') 3166 end if 3167 if(Sigp%nomegasr==1) then 3168 domegas=0.d0 3169 else 3170 domegas=(Sigp%maxomega_r-Sigp%minomega_r)/(Sigp%nomegasr-1) 3171 endif 3172 !TODO this should be moved to Sr% and done in init_sigma_t 3173 ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr)) 3174 do io=1,Sigp%nomegasr 3175 Sigp%omega_r(io) = CMPLX(Sigp%minomega_r + domegas*(io-1),zero) 3176 end do 3177 write(msg,'(4a,i8,3(2a,f8.3),a)')ch10,& 3178 ' Parameters for the calculation of the spectral function : ',ch10,& 3179 ' Number of points = ',Sigp%nomegasr,ch10,& 3180 ' Min frequency [eV] = ',Sigp%minomega_r*Ha_eV,ch10,& 3181 ' Max frequency [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,& 3182 ' Frequency step [eV] = ',domegas*Ha_eV,ch10 3183 call wrtout(std_out, msg) 3184 else 3185 Sigp%minomega_r = MINVAL(Dtset%gw_freqsp(:)) 3186 Sigp%maxomega_r = MAXVAL(Dtset%gw_freqsp(:)) 3187 !TODO this should be moved to Sr% and done in init_sigma_t 3188 ABI_MALLOC(Sigp%omega_r,(Sigp%nomegasr)) 3189 do io=1,Sigp%nomegasr 3190 Sigp%omega_r(io) = CMPLX(Dtset%gw_freqsp(io),zero) 3191 end do 3192 write(msg,'(4a,i8,2(2a,f8.3),3a)')ch10,& 3193 ' Parameters for the calculation of the spectral function : ',ch10,& 3194 ' Number of points = ',Sigp%nomegasr,ch10,& 3195 ' Min frequency [eV] = ',Sigp%minomega_r*Ha_eV,ch10,& 3196 ' Max frequency [eV] = ',Sigp%maxomega_r*Ha_eV,ch10,& 3197 ' A custom set of frequencies is used! See the input file for values.',ch10 3198 call wrtout(std_out, msg) 3199 end if 3200 else 3201 !In indefo all these quantities are set to zero 3202 !Sigp%nomegasr=1 3203 !allocate(Sigp%omega_r(Sigp%nomegasr)) 3204 !Sigp%omega_r(1)=0 3205 end if 3206 3207 ! Dimensional primitive translations rprimd (from input), gprimd, metrics and unit cell volume 3208 call mkrdim(acell,rprim,rprimd) 3209 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 3210 3211 Sigp%npwwfn = Dtset%npwwfn 3212 Sigp%npwx = Dtset%npwsigx 3213 3214 ! Read parameters of the WFK, verifify them and retrieve all G-vectors. 3215 call wfk_read_eigenvalues(wfk_fname,energies_p,Hdr_wfk,comm) 3216 mband = MAXVAL(Hdr_wfk%nband) 3217 3218 remove_inv = .FALSE. 3219 call hdr_wfk%vs_dtset(dtset) 3220 3221 test_npwkss = 0 3222 call make_gvec_kss(Dtset%nkpt,Dtset%kptns,Hdr_wfk%ecut_eff,Dtset%symmorphi,Dtset%nsym,Dtset%symrel,Dtset%tnons,& 3223 gprimd,Dtset%prtvol,test_npwkss,test_gvec_kss,ierr) 3224 ABI_CHECK(ierr==0,"Fatal error in make_gvec_kss") 3225 3226 ABI_MALLOC(gvec_kss,(3, test_npwkss)) 3227 gvec_kss = test_gvec_kss 3228 ng_kss = test_npwkss 3229 3230 ng = MIN(SIZE(gvec_kss,DIM=2),SIZE(test_gvec_kss,DIM=2)) 3231 ierr = 0 3232 do ig1=1,ng 3233 if (ANY(gvec_kss(:,ig1)/=test_gvec_kss(:,ig1))) then 3234 ierr=ierr+1 3235 write(std_out,*)" gvec_kss ",ig1,"/",ng,gvec_kss(:,ig1),test_gvec_kss(:,ig1) 3236 end if 3237 end do 3238 ABI_CHECK(ierr == 0, "Mismatch between gvec_kss and test_gvec_kss") 3239 3240 ABI_FREE(test_gvec_kss) 3241 3242 ! Get important dimensions from the WFK header 3243 Sigp%nsppol = Hdr_wfk%nsppol 3244 Sigp%nspinor = Hdr_wfk%nspinor 3245 Sigp%nsig_ab = Hdr_wfk%nspinor**2 ! TODO Is it useful calculating only diagonal terms? 3246 3247 if (Sigp%nbnds > mband) then 3248 Sigp%nbnds = mband 3249 Dtset%nband(:) = mband 3250 Dtset%mband = MAXVAL(Dtset%nband) 3251 write(msg,'(3a,i4,a)')& 3252 'Number of bands found less then required',ch10,& 3253 'calculation will proceed with nbnds = ',mband,ch10 3254 ABI_WARNING(msg) 3255 end if 3256 3257 ! Check input 3258 if (Sigp%ppmodel==3.or.Sigp%ppmodel==4) then 3259 if (gwcalctyp>=10) then 3260 write(msg,'(a,i3,a)')' The ppmodel chosen and gwcalctyp ',Dtset%gwcalctyp,' are not compatible. ' 3261 ABI_ERROR(msg) 3262 end if 3263 if (Sigp%nspinor==2) then 3264 write(msg,'(a,i3,a)')' The ppmodel chosen and nspinor ',Sigp%nspinor,' are not compatible. ' 3265 ABI_ERROR(msg) 3266 end if 3267 end if 3268 3269 ! Create crystal_t data type 3270 cryst = Hdr_wfk%get_crystal(gw_timrev=timrev, remove_inv=remove_inv) 3271 call cryst%print() 3272 3273 if (Sigp%npwwfn > ng_kss) then ! cannot use more G"s for the wfs than those stored on file 3274 Sigp%npwwfn =ng_kss 3275 Dtset%npwwfn =ng_kss 3276 write(msg,'(2a,(a,i8,a))')& 3277 'Number of G-vectors for WFS found in the KSS file is less than required',ch10,& 3278 'calculation will proceed with npwwfn = ',Sigp%npwwfn,ch10 3279 ABI_WARNING(msg) 3280 end if 3281 3282 if (Sigp%npwx>ng_kss) then 3283 ! Have to recalcuate the (large) sphere for Sigma_x. 3284 pinv=1; if (remove_inv.and.Cryst%timrev==2) pinv=-1 3285 gamma_point(:,1) = (/zero,zero,zero/); nullify(gsphere_sigx_p) 3286 3287 call merge_and_sort_kg(1,gamma_point,Dtset%ecutsigx,Cryst%nsym,pinv,Cryst%symrel,& 3288 Cryst%gprimd,gsphere_sigx_p,Dtset%prtvol) 3289 3290 ng_sigx=SIZE(gsphere_sigx_p,DIM=2) 3291 Sigp%npwx = ng_sigx 3292 Dtset%npwsigx = ng_sigx 3293 3294 write(msg,'(2a,(a,i8,a))')& 3295 'Number of G-vectors for Sigma_x found in the KSS file is less than required',ch10,& 3296 'calculation will proceed with npwsigx = ',Sigp%npwx,ch10 3297 ABI_WARNING(msg) 3298 3299 ltest = (Sigp%npwx >= ng_kss) 3300 ABI_CHECK(ltest,"Sigp%npwx<ng_kss!") 3301 3302 ! * Fill gvec_kss with larger sphere. 3303 ABI_FREE(gvec_kss) 3304 ABI_MALLOC(gvec_kss,(3,Sigp%npwx)) 3305 gvec_kss = gsphere_sigx_p 3306 ABI_FREE(gsphere_sigx_p) 3307 end if 3308 3309 ! Set up of the k-points and tables in the whole BZ === 3310 ! TODO Recheck symmorphy and inversion 3311 call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.FALSE.) 3312 !call Kmesh%init(Cryst,Hdr_wfk%nkpt,Hdr_wfk%kptns,Dtset%kptopt,wrap_1zone=.TRUE.) 3313 3314 ! Some required information are not filled up inside kmesh_init 3315 ! So doing it here, even though it is not clean 3316 Kmesh%kptrlatt(:,:) =Dtset%kptrlatt(:,:) 3317 Kmesh%nshift =Dtset%nshiftk 3318 ABI_MALLOC(Kmesh%shift,(3,Kmesh%nshift)) 3319 Kmesh%shift(:,:) =Dtset%shiftk(:,1:Dtset%nshiftk) 3320 3321 call Kmesh%print("K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL") 3322 call Kmesh%print("K-mesh for the wavefunctions",ab_out, 0, "COLL") 3323 3324 ! === Initialize the band structure datatype === 3325 ! * Copy WFK energies and occupations up to Sigp%nbnds==Dtset%nband(:) 3326 bantot = SUM(Dtset%nband(1:Dtset%nkpt*nsppol)) 3327 ABI_MALLOC(doccde,(bantot)) 3328 ABI_MALLOC(eigen,(bantot)) 3329 ABI_MALLOC(occfact,(bantot)) 3330 doccde(:)=zero; eigen(:)=zero; occfact(:)=zero 3331 3332 jj=0; ibtot=0 3333 do isppol=1,nsppol 3334 do ikibz=1,Dtset%nkpt 3335 do ib=1,Hdr_wfk%nband(ikibz+(isppol-1)*Dtset%nkpt) 3336 ibtot=ibtot+1 3337 if (ib<=Sigp%nbnds) then 3338 jj=jj+1 3339 occfact(jj)=Hdr_wfk%occ(ibtot) 3340 eigen (jj)=energies_p(ib,ikibz,isppol) 3341 end if 3342 end do 3343 end do 3344 end do 3345 ABI_FREE(energies_p) 3346 3347 ! Make sure that Dtset%wtk==Kmesh%wt due to the dirty treatment of 3348 ! symmetry operations in the old GW code (symmorphy and inversion) 3349 ltest=(ALL(ABS(Dtset%wtk(1:Kmesh%nibz)-Kmesh%wt(1:Kmesh%nibz))<tol6)) 3350 ABI_CHECK(ltest,'Mismatch between Dtset%wtk and Kmesh%wt') 3351 3352 ABI_MALLOC(npwarr,(Dtset%nkpt)) 3353 npwarr(:)=Sigp%npwwfn 3354 3355 call ebands_init(bantot,ks_ebands,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,& 3356 doccde,eigen,Dtset%istwfk,Kmesh%ibz,Dtset%nband,& 3357 Kmesh%nibz,npwarr,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt,& 3358 dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig,& 3359 dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 3360 3361 ABI_FREE(doccde) 3362 ABI_FREE(eigen) 3363 ABI_FREE(npwarr) 3364 3365 ! Calculate KS occupation numbers and ks_vbk(nkibz, nsppol) ==== 3366 ! ks_vbk gives the (valence|last Fermi band) index for each k and spin. 3367 ! spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals 3368 call ebands_update_occ(ks_ebands, Dtset%spinmagntarget, prtvol=0) 3369 3370 gaps = ebands_get_gaps(ks_ebands, gap_err) 3371 call gaps%print(unit=std_out) 3372 call ebands_report_gap(ks_ebands, unit=std_out) 3373 3374 ABI_MALLOC(val_indices,(ks_ebands%nkpt, nsppol)) 3375 val_indices = ebands_get_valence_idx(ks_ebands) 3376 3377 ! Create Sigma header 3378 ! TODO Fix problems with symmorphy and k-points 3379 call hdr_init(ks_ebands,codvsn,Dtset,Hdr_out,Pawtab,pertcase0,Psps,wvl) 3380 3381 ! Get Pawrhoij from the header of the WFK file 3382 ABI_MALLOC(Pawrhoij, (Cryst%natom*Dtset%usepaw)) 3383 if (Dtset%usepaw==1) then 3384 call pawrhoij_alloc(Pawrhoij, 1, Dtset%nspden, Dtset%nspinor, nsppol, Cryst%typat, pawtab=Pawtab) 3385 call pawrhoij_copy(Hdr_wfk%Pawrhoij, Pawrhoij) 3386 end if 3387 3388 call hdr_out%update(bantot,1.0d20,1.0d20,1.0d20,1.0d20,Cryst%rprimd,occfact,Pawrhoij,Cryst%xred,dtset%amu_orig(:,1)) 3389 ABI_FREE(occfact) 3390 call pawrhoij_free(Pawrhoij) 3391 ABI_FREE(Pawrhoij) 3392 3393 ! =========================================================== 3394 ! ==== Setup of k-points and bands for the GW corrections ==== 3395 ! =========================================================== 3396 ! * maxbdgw and minbdgw are the Max and min band index for GW corrections over k-points. 3397 ! They are used to dimension the wavefunctions and to calculate the matrix elements. 3398 ! 3399 if (dtset%nkptgw == 0) then 3400 ! 3401 ! Use qp_range to select the interesting k-points and the corresponing bands. 3402 ! 3403 ! 0 --> Compute the QP corrections only for the fundamental and the direct gap. 3404 ! +num --> Compute the QP corrections for all the k-points in the irreducible zone and include `num` 3405 ! bands above and below the Fermi level. 3406 ! -num --> Compute the QP corrections for all the k-points in the irreducible zone. 3407 ! Include all occupied states and `num` empty states. 3408 3409 call wrtout(std_out, "nkptgw == 0 ==> Automatic selection of k-points and bands for the corrections.") 3410 3411 if (gap_err /= 0 .and. dtset%gw_qprange == 0) then 3412 msg = "Problem while computing the fundamental and direct gap (likely metal). Will replace gw_qprange=0 with gw_qprange=1" 3413 ABI_WARNING(msg) 3414 dtset%gw_qprange = 1 3415 end if 3416 3417 gw_qprange = dtset%gw_qprange 3418 3419 if (dtset%ucrpa > 0) then 3420 dtset%nkptgw = Kmesh%nbz 3421 Sigp%nkptgw = dtset%nkptgw 3422 ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw)) 3423 ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol)) 3424 ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol)) 3425 Sigp%kptgw(:,:) = Kmesh%bz(:,:) 3426 Sigp%minbnd = 1 3427 Sigp%maxbnd = Sigp%nbnds 3428 3429 else if (gw_qprange /= 0) then 3430 ! Include all the k-points in the IBZ. 3431 ! Note that kptgw == ebands%kptns so we can use a single ik index in the loop over k-points 3432 ! No need to map kptgw onto ebands%kptns. 3433 dtset%nkptgw = Kmesh%nibz 3434 Sigp%nkptgw = dtset%nkptgw 3435 ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw)) 3436 ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol)) 3437 ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol)) 3438 Sigp%kptgw(:,:) = Kmesh%ibz(:,:) 3439 Sigp%minbnd = 1 3440 Sigp%maxbnd = Sigp%nbnds 3441 3442 if (gw_qprange > 0) then 3443 ! All k-points: Add buffer of bands above and below the Fermi level. 3444 do spin=1,nsppol 3445 do ik=1,Sigp%nkptgw 3446 Sigp%minbnd(ik, spin) = MAX(val_indices(ik,spin) - gw_qprange, 1) 3447 Sigp%maxbnd(ik, spin) = MIN(val_indices(ik,spin) + gw_qprange + 1, Sigp%nbnds) 3448 end do 3449 end do 3450 3451 else 3452 ! All k-points: include all occupied states and -gw_qprange empty states. 3453 Sigp%minbnd = 1 3454 do spin=1,nsppol 3455 do ik=1,Sigp%nkptgw 3456 Sigp%maxbnd(ik, spin) = MIN(val_indices(ik,spin) - gw_qprange, Sigp%nbnds) 3457 end do 3458 end do 3459 end if 3460 3461 else 3462 ! gw_qprange is not specified in the input. 3463 ! Include the direct and the fundamental KS gap. 3464 ! The main problem here is that kptgw and nkptgw do not depend on the spin and therefore 3465 ! we have compute the union of the k-points where the fundamental and the direct gaps are located. 3466 ! 3467 ! Find the list of `interesting` kpoints. 3468 ABI_CHECK(gap_err == 0, "gw_qprange 0 cannot be used because I cannot find the gap (gap_err !=0)") 3469 nk_found = 1; kpos(1) = gaps%fo_kpos(1,1) 3470 3471 do spin=1,nsppol 3472 do ifo=1,3 3473 ik = gaps%fo_kpos(ifo, spin) 3474 found = .FALSE.; jj = 0 3475 do while (.not. found .and. jj < nk_found) 3476 jj = jj + 1; found = (kpos(jj) == ik) 3477 end do 3478 if (.not. found) then 3479 nk_found = nk_found + 1; kpos(nk_found) = ik 3480 end if 3481 end do 3482 end do 3483 3484 ! Now we can define the list of k-points and the bands range. 3485 dtset%nkptgw = nk_found 3486 Sigp%nkptgw = dtset%nkptgw 3487 3488 ABI_MALLOC(Sigp%kptgw,(3, Sigp%nkptgw)) 3489 ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol)) 3490 ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol)) 3491 3492 do ii=1,Sigp%nkptgw 3493 ik = kpos(ii) 3494 Sigp%kptgw(:,ii) = Kmesh%ibz(:,ik) 3495 do spin=1,nsppol 3496 Sigp%minbnd(ii,spin) = val_indices(ik, spin) 3497 Sigp%maxbnd(ii,spin) = val_indices(ik, spin) + 1 3498 end do 3499 end do 3500 end if 3501 3502 else 3503 ! Treat only the k-points and bands specified in the input file. 3504 Sigp%nkptgw = dtset%nkptgw 3505 ABI_MALLOC(Sigp%kptgw, (3, Sigp%nkptgw)) 3506 ABI_MALLOC(Sigp%minbnd, (Sigp%nkptgw, nsppol)) 3507 ABI_MALLOC(Sigp%maxbnd, (Sigp%nkptgw, nsppol)) 3508 3509 do spin=1,nsppol 3510 Sigp%minbnd(:,spin)= dtset%bdgw(1,:,spin) 3511 Sigp%maxbnd(:,spin)= dtset%bdgw(2,:,spin) 3512 end do 3513 3514 do ii=1,3 3515 do ikcalc=1,Sigp%nkptgw 3516 Sigp%kptgw(ii,ikcalc) = Dtset%kptgw(ii,ikcalc) 3517 end do 3518 end do 3519 3520 do spin=1,nsppol 3521 do ikcalc=1,Sigp%nkptgw 3522 if (Dtset%bdgw(2,ikcalc,spin) > Sigp%nbnds) then 3523 write(msg,'(a,2i0,2(a,i0),2a,i0)')& 3524 "For (k,s) ",ikcalc,spin," bdgw= ",Dtset%bdgw(2,ikcalc,spin), " > nbnds=",Sigp%nbnds,ch10,& 3525 "Calculation will continue with bdgw =",Sigp%nbnds 3526 ABI_COMMENT(msg) 3527 Dtset%bdgw(2, ikcalc, spin) = Sigp%nbnds 3528 end if 3529 end do 3530 end do 3531 3532 end if 3533 3534 ! Make sure that all the degenerate states are included. 3535 ! * We will have to average the GW corrections over degenerate states if symsigma=1 is used. 3536 ! * KS states belonging to the same irreducible representation should be included in the basis set used for SCGW. 3537 if (Sigp%symsigma /=0 .or. gwcalctyp >= 10) then 3538 do isppol=1,nsppol 3539 do ikcalc=1,Sigp%nkptgw 3540 3541 if (kmesh%has_IBZ_item(Sigp%kptgw(:,ikcalc), ikibz, G0)) then 3542 call ebands_enclose_degbands(ks_ebands,ikibz,isppol, & 3543 Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol),changed,tol_enediff) 3544 3545 if (changed) then 3546 write(msg,'(2(a,i0),2a,2(1x,i0))')& 3547 "Not all the degenerate states at ikcalc= ",ikcalc,", spin= ",isppol, ch10, & 3548 "were included in the bdgw set. bdgw has been changed to: ",Sigp%minbnd(ikcalc,isppol),Sigp%maxbnd(ikcalc,isppol) 3549 ABI_COMMENT(msg) 3550 end if 3551 else 3552 write(msg,'(3a)')& 3553 ' not in the list of k points treated in the preparatory SCF run.',ch10, & 3554 ' Change kptgw, or shiftk of previous run.' 3555 ABI_ERROR(sjoin('k-point', ktoa(Sigp%kptgw(:,ikcalc)),trim(msg))) 3556 end if 3557 3558 end do 3559 end do 3560 end if 3561 3562 Sigp%minbdgw = MINVAL(Sigp%minbnd) 3563 Sigp%maxbdgw = MAXVAL(Sigp%maxbnd) 3564 3565 ! Check if there are duplicated k-point in Sigp% 3566 do ii=1,Sigp%nkptgw 3567 do jj=ii+1,Sigp%nkptgw 3568 if (isamek(Sigp%kptgw(:,ii), Sigp%kptgw(:,jj), G0)) then 3569 write(msg,'(5a)')& 3570 'kptgw contains duplicated k-points. This is not allowed since ',ch10,& 3571 'the QP corrections for this k-point will be calculated more than once. ',ch10,& 3572 'Check your input file. ' 3573 ABI_ERROR(msg) 3574 end if 3575 end do 3576 end do 3577 3578 !=== Check if the k-points are in the BZ === 3579 ! FB: Honestly the code is not able to treat k-points, which are not in the input list of points in the IBZ. 3580 ! This extension should require to change the code in different places. 3581 ! Therefore, one should by now prevent the user from calculating sigma for a k-point not in the IBZ. 3582 ABI_MALLOC(Sigp%kptgw2bz, (Sigp%nkptgw)) 3583 3584 do ikcalc=1,Sigp%nkptgw 3585 if (kmesh%has_BZ_item(Sigp%kptgw(:,ikcalc), ikcalc2bz, G0)) then 3586 Sigp%kptgw2bz(ikcalc) = ikcalc2bz 3587 else 3588 write(msg,'(3a)')& 3589 ' not in the list of k points treated in the preparatory SCF run.',ch10,' Change kptgw, or shiftk of previous run.' 3590 ABI_ERROR(sjoin('k-point:', ktoa(Sigp%kptgw(:,ikcalc)),trim(msg))) 3591 end if 3592 end do 3593 3594 ! Warn the user if SCGW run and not all the k-points are included. 3595 if (gwcalctyp >= 10 .and. Sigp%nkptgw /= Hdr_wfk%nkpt) then 3596 write(msg,'(3a,2(a,i0),2a)')ch10,& 3597 " COMMENT: In a self-consistent GW run, the QP corrections should be calculated for all the k-points of the KSS file ",ch10,& 3598 " but nkptgw= ",Sigp%nkptgw," and WFK nkpt= ",Hdr_wfk%nkpt,ch10,& 3599 " Assuming expert user. Execution will continue. " 3600 call wrtout(ab_out, msg) 3601 end if 3602 3603 ! Setup of the table used in the case of SCGW on wavefunctions to reduce the number 3604 ! of elements <i,kgw,s|\Sigma|j,kgw,s> that have to be calculated. No use of symmetries, except for Hermiticity. 3605 call sigma_tables(Sigp, Kmesh) 3606 3607 ! === Read external file and initialize basic dimension of Er% === 3608 ! TODO use mqmem as input variable instead of gwmem 3609 3610 ! === If required, use a matrix for $\Sigma_c$ which is smaller than that stored on file === 3611 ! * By default the entire matrix is read and used, 3612 ! * Define consistently npweps and ecuteps for \Sigma_c according the input 3613 if (Dtset%npweps>0.or.Dtset%ecuteps>0) then 3614 if (Dtset%npweps>0) Dtset%ecuteps=zero 3615 nsheps=0 3616 call setshells(Dtset%ecuteps,Dtset%npweps,nsheps,Dtset%nsym,gmet,gprimd,Dtset%symrel,'eps',ucvol) 3617 end if 3618 3619 mqmem=0; if (Dtset%gwmem/10==1) mqmem=1 3620 3621 if (dtset%getscr /=0 .or. dtset%irdscr/=0 .or. dtset%getscr_filepath /= ABI_NOFILE) then 3622 fname=Dtfil%fnameabi_scr 3623 else if (Dtset%getsuscep/=0.or.Dtset%irdsuscep/=0) then 3624 fname=Dtfil%fnameabi_sus 3625 else 3626 fname=Dtfil%fnameabi_scr 3627 !FIXME this has to be cleaned, in tgw2_3 Dtset%get* and Dtset%ird* are not defined 3628 !ABI_ERROR("getsuscep or irdsuscep are not defined") 3629 end if 3630 ! 3631 ! === Setup of q-mesh in the whole BZ === 3632 ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 3633 ! epsilon^-1(Sq) should be symmetrized in csigme using a different expression (G-G_o is needed) 3634 ! 3635 if (sigma_needs_w(Sigp)) then 3636 if (.not. file_exists(fname)) then 3637 fname = nctk_ncify(fname) 3638 ABI_COMMENT(sjoin("File not found. Will try netcdf file:", fname)) 3639 end if 3640 3641 call init_Er_from_file(Er,fname,mqmem,Dtset%npweps,comm) 3642 3643 Sigp%npwc=Er%npwe 3644 if (Sigp%npwc>Sigp%npwx) then 3645 Sigp%npwc=Sigp%npwx 3646 ABI_COMMENT("Found npw_correlation > npw_exchange, Imposing npwc=npwx") 3647 ! There is a good reason for doing so, see csigme.F90 and the size of the arrays 3648 ! rhotwgp and rhotwgp: we need to define a max size and we opt for Sigp%npwx. 3649 end if 3650 Er%npwe=Sigp%npwc 3651 Dtset%npweps=Er%npwe 3652 call Qmesh%init(Cryst,Er%nqibz,Er%qibz,Dtset%kptopt) 3653 3654 else 3655 Er%npwe =1 3656 Sigp%npwc =1 3657 Dtset%npweps=1 3658 call find_qmesh(Qmesh,Cryst,Kmesh) 3659 ABI_MALLOC(Er%gvec,(3,1)) 3660 Er%gvec(:,1) = [0, 0, 0] 3661 end if 3662 3663 call Qmesh%print("Q-mesh for screening function",std_out,Dtset%prtvol,"COLL") 3664 call Qmesh%print("Q-mesh for screening function",ab_out ,0 ,"COLL") 3665 3666 3667 do iqbz=1,Qmesh%nbz 3668 call qmesh%get_BZ_item(iqbz, q_bz, iq_ibz, isym, itim, umklp=q_umklp) 3669 !print *, "iqbz, q_bz, iq_ibz, isym, itim, umklp" 3670 !print *, iqbz, q_bz, iq_ibz, isym, itim, q_umklp 3671 3672 if (ANY(q_umklp/=0)) then 3673 sq = (3-2*itim)*MATMUL(Cryst%symrec(:,:,isym),Qmesh%ibz(:,iq_ibz)) 3674 write(std_out,*) sq,Qmesh%bz(:,iqbz) 3675 write(msg,'(a,3f6.3,a,3f6.3,2a,9i3,a,i2,2a)')& 3676 'qpoint ',Qmesh%bz(:,iqbz),' is the symmetric of ',Qmesh%ibz(:,iq_ibz),ch10,& 3677 'through operation ',Cryst%symrec(:,:,isym),' and itim ',itim,ch10,& 3678 'however a non-zero umklapp G_o vector is required and this is not yet allowed' 3679 ABI_ERROR(msg) 3680 end if 3681 end do 3682 !stop 3683 ! 3684 ! === Find optimal value for G-sphere enlargment due to oscillator matrix elements === 3685 ! * Here I have to be sure that Qmesh%bz is always inside the BZ, not always true size bz is buggy 3686 ! * -one is used because we loop over all the possibile differences, unlike screening 3687 3688 call get_ng0sh(Sigp%nkptgw,Sigp%kptgw,Kmesh%nbz,Kmesh%bz,Qmesh%nbz,Qmesh%bz,-one,ng0sh_opt) 3689 call wrtout(std_out, sjoin(' Optimal value for ng0sh ', ltoa(ng0sh_opt))) 3690 Sigp%mG0 = ng0sh_opt 3691 3692 ! G-sphere for W and Sigma_c is initialized from the SCR file. 3693 call Gsph_c%init(Cryst, Er%npwe, gvec=Er%gvec) 3694 call Gsph_x%init(Cryst, Sigp%npwx, gvec=gvec_kss) 3695 Sigp%ecuteps = Gsph_c%ecut 3696 Dtset%ecuteps = Sigp%ecuteps 3697 3698 ! === Make biggest G-sphere of Sigp%npwvec vectors === 3699 Sigp%npwvec=MAX(Sigp%npwwfn,Sigp%npwx) 3700 call Gsph_Max%init(Cryst, Sigp%npwvec, gvec=gvec_kss) 3701 !BEGINDEBUG 3702 ! Make sure that the two G-spheres are equivalent. 3703 ierr=0 3704 if (sigma_needs_w(Sigp)) then 3705 ng = MIN(SIZE(Gsph_c%gvec,DIM=2),SIZE(gvec_kss,DIM=2)) 3706 do ig1=1,ng 3707 if (ANY(Gsph_c%gvec(:,ig1)/=gvec_kss(:,ig1))) then 3708 ierr=ierr+1 3709 write(std_out,*)" Gsph_c, gvec_kss ",ig1,"/",ng,Gsph_c%gvec(:,ig1),gvec_kss(:,ig1) 3710 end if 3711 end do 3712 ABI_CHECK(ierr==0,"Mismatch between Gsph_c and gvec_kss") 3713 end if 3714 ierr=0 3715 ng = MIN(SIZE(Gsph_x%gvec,DIM=2),SIZE(gvec_kss,DIM=2)) 3716 do ig1=1,ng 3717 if (ANY(Gsph_x%gvec(:,ig1)/=gvec_kss(:,ig1))) then 3718 ierr=ierr+1 3719 write(std_out,*)" Gsph_x, gvec_kss ",ig1,"/",ng,Gsph_x%gvec(:,ig1),gvec_kss(:,ig1) 3720 end if 3721 end do 3722 ABI_CHECK(ierr==0,"Mismatch between Gsph_x and gvec_kss") 3723 !ENDDEBUG 3724 3725 ABI_FREE(gvec_kss) 3726 ! 3727 ! === Get Fourier components of the Coulomb term for all q-points in the IBZ === 3728 ! * If required, use a cutoff in the interaction 3729 ! * Pcv%vc_sqrt contains Vc^{-1/2} 3730 ! * Setup also the analytical calculation of the q->0 component 3731 ! FIXME recheck ngfftf since I got different charge outside the cutoff region 3732 3733 if (Dtset%gw_nqlwl==0) then 3734 nqlwl=1 3735 ABI_MALLOC(qlwl,(3,nqlwl)) 3736 qlwl(:,1)= GW_Q0_DEFAULT 3737 else 3738 nqlwl=Dtset%gw_nqlwl 3739 ABI_MALLOC(qlwl,(3,nqlwl)) 3740 qlwl(:,:)=Dtset%gw_qlwl(:,1:nqlwl) 3741 end if 3742 3743 ! The Coulomb interaction used here might have two terms: 3744 ! the first term generates the usual sigma self-energy, but possibly, one should subtract 3745 ! from it the Coulomb interaction already present in the Kohn-Sham basis, 3746 ! if the usefock associated to ixc is one. 3747 ! The latter excludes (in the present implementation) mod(Dtset%gwcalctyp,10)==5 3748 nvcoul_init=1 3749 call get_xclevel(Dtset%ixc,xclevel_ixc,usefock_ixc) 3750 3751 if(usefock_ixc==1)then 3752 nvcoul_init=2 3753 if(mod(Dtset%gwcalctyp,10)==5)then 3754 write(msg,'(4a,i3,a,i3,4a,i5)')ch10,& 3755 ' The starting wavefunctions were obtained from self-consistent calculations in the planewave basis set',ch10,& 3756 ' with ixc = ',Dtset%ixc,' associated with usefock =',usefock_ixc,ch10,& 3757 ' In this case, the present implementation does not allow that the self-energy for sigma corresponds to',ch10,& 3758 ' mod(gwcalctyp,10)==5, while your gwcalctyp= ',Dtset%gwcalctyp 3759 ABI_ERROR(msg) 3760 endif 3761 endif 3762 do ivcoul_init=1,nvcoul_init 3763 rcut = Dtset%rcut 3764 icutcoul_eff=Dtset%gw_icutcoul 3765 Sigp%sigma_mixing=one 3766 if( mod(Dtset%gwcalctyp,10)==5 .or. ivcoul_init==2)then 3767 if(abs(Dtset%hyb_mixing)>tol8)then 3768 ! Warning: the absolute value is needed, because of the singular way used to define the default for this input variable 3769 Sigp%sigma_mixing=abs(Dtset%hyb_mixing) 3770 else if(abs(Dtset%hyb_mixing_sr)>tol8)then 3771 Sigp%sigma_mixing=abs(Dtset%hyb_mixing_sr) 3772 icutcoul_eff=5 3773 endif 3774 if(abs(rcut)<tol6 .and. abs(Dtset%hyb_range_fock)>tol8)rcut=one/Dtset%hyb_range_fock 3775 endif 3776 3777 if (ivcoul_init == 1) then 3778 3779 if (Gsph_x%ng > Gsph_c%ng) then 3780 call Vcp%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 3781 Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm) 3782 else 3783 call Vcp%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 3784 Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm) 3785 end if 3786 3787 else 3788 3789 ! Use a temporary Vcp_ks to compute the Coulomb interaction already present 3790 ! in the Fock part of the Kohn-Sham Hamiltonian 3791 if (Gsph_x%ng > Gsph_c%ng) then 3792 call Vcp_ks%init(Gsph_x,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 3793 Dtset%ecutsigx,Gsph_x%ng,nqlwl,qlwl,comm) 3794 else 3795 call Vcp_ks%init(Gsph_c,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,Dtset%vcutgeo,& 3796 Dtset%ecutsigx,Gsph_c%ng,nqlwl,qlwl,comm) 3797 end if 3798 3799 ! Now compute the residual Coulomb interaction. 3800 Vcp%vc_sqrt_resid=sqrt(Vcp%vc_sqrt**2-Sigp%sigma_mixing*Vcp_ks%vc_sqrt**2) 3801 Vcp%i_sz_resid=Vcp%i_sz-Sigp%sigma_mixing*Vcp_ks%i_sz 3802 3803 ! The mixing factor has already been accounted for, so set it back to one 3804 Sigp%sigma_mixing=one 3805 call Vcp_ks%free() 3806 3807 !write(std_out,'(a)')' setup_sigma : the residual Coulomb interaction has been computed' 3808 endif 3809 3810 !#else 3811 ! call Vcp%init(Gsph_Max,Cryst,Qmesh,Kmesh,rcut,icutcoul_eff,ivcoul_init,Dtset%vcutgeo,& 3812 !& Dtset%ecutsigx,Sigp%npwx,nqlwl,qlwl,ngfftf,comm) 3813 !#endif 3814 3815 end do 3816 3817 #if 0 3818 ! Using the random q for the optical limit is one of the reasons 3819 ! why sigma breaks the initial energy degeneracies. 3820 Vcp%i_sz=zero 3821 Vcp%vc_sqrt(1,:)=czero 3822 Vcp%vcqlwl_sqrt(1,:)=czero 3823 #endif 3824 3825 ABI_FREE(qlwl) 3826 3827 Sigp%ecuteps = Dtset%ecuteps 3828 Sigp%ecutwfn = Dtset%ecutwfn 3829 Sigp%ecutsigx = Dtset%ecutsigx 3830 3831 ! === Setup of the FFT mesh for the oscilator strengths === 3832 ! * Init gwc_ngfft(7:18) and gwx_ngfft(7:18) with Dtset%ngfft(7:18) 3833 ! * Here we redefine gwc_ngfft(1:6) according to the following options: 3834 ! 3835 ! method == 0 --> FFT grid read from fft.in (debugging purpose) 3836 ! method == 1 --> Normal FFT mesh 3837 ! method == 2 --> Slightly augmented FFT grid to calculate exactly rho_tw_g (see setmesh.F90) 3838 ! method == 3 --> Doubled FFT grid, same as the the FFT for the density, 3839 ! 3840 ! enforce_sym == 1 --> Enforce a FFT mesh compatible with all the symmetry operation and FFT library 3841 ! enforce_sym == 0 --> Find the smallest FFT grid compatbile with the library, do not care about symmetries 3842 ! 3843 gwc_ngfft(1:18) = Dtset%ngfft(1:18) 3844 gwx_ngfft(1:18) = Dtset%ngfft(1:18) 3845 3846 method = 2 3847 if (Dtset%fftgw == 00 .or. Dtset%fftgw == 01) method = 0 3848 if (Dtset%fftgw == 10 .or. Dtset%fftgw == 11) method = 1 3849 if (Dtset%fftgw == 20 .or. Dtset%fftgw == 21) method = 2 3850 if (Dtset%fftgw == 30 .or. Dtset%fftgw == 31) method = 3 3851 enforce_sym = mod(dtset%fftgw, 10) 3852 3853 ! FFT mesh for sigma_x. 3854 call setmesh(gmet, Gsph_Max%gvec, gwx_ngfft, Sigp%npwvec, Sigp%npwx, Sigp%npwwfn, & 3855 gwx_nfftot, method, Sigp%mG0, Cryst, enforce_sym) 3856 3857 ! FFT mesh for sigma_c. 3858 call setmesh(gmet, Gsph_Max%gvec, gwc_ngfft, Sigp%npwvec, Er%npwe, Sigp%npwwfn,& 3859 gwc_nfftot, method, Sigp%mG0, Cryst, enforce_sym, unit=dev_null) 3860 3861 ! ====================================================================== 3862 ! ==== Check for presence of files with core orbitals, for PAW only ==== 3863 ! ====================================================================== 3864 Sigp%use_sigxcore=0 3865 if (Dtset%usepaw==1.and.Dtset%gw_sigxcore==1) then 3866 ii = 0 3867 do itypat=1,Cryst%ntypat 3868 string = Psps%filpsp(itypat) 3869 fcore = "CORE_"//TRIM(basename(string)) 3870 ic = INDEX (TRIM(string), "/" , back=.TRUE.) 3871 if (ic>0 .and. ic<LEN_TRIM(string)) then 3872 ! string defines a path, prepend path to fcore 3873 fcore = Psps%filpsp(itypat)(1:ic)//TRIM(fcore) 3874 end if 3875 if (file_exists(fcore)) then 3876 ii = ii+1 3877 else 3878 ABI_WARNING(sjoin("HF decoupling is required but cannot find file:", fcore)) 3879 end if 3880 end do 3881 3882 Sigp%use_sigxcore=1 3883 if (ii/=Cryst%ntypat) then 3884 ABI_ERROR("Files with core orbitals not found") 3885 end if 3886 end if ! PAW+HF decoupling 3887 ! 3888 ! ============================== 3889 ! ==== Extrapolar technique ==== 3890 ! ============================== 3891 Sigp%gwcomp = Dtset%gwcomp 3892 Sigp%gwencomp = Dtset%gwencomp 3893 3894 if (Sigp%gwcomp==1) then 3895 write(msg,'(6a,e11.4,a)')ch10,& 3896 'Using the extrapolar approximation to accelerate convergence',ch10,& 3897 'with respect to the number of bands included',ch10,& 3898 'with gwencomp: ',Sigp%gwencomp*Ha_eV,' [eV]' 3899 call wrtout(std_out, msg) 3900 end if 3901 ! 3902 ! =================================== 3903 ! ==== Final compatibility tests ==== 3904 ! =================================== 3905 ltest=(ks_ebands%mband == Sigp%nbnds .and. ALL(ks_ebands%nband == Sigp%nbnds)) 3906 ABI_CHECK(ltest,'BUG in definition of ks_ebands%nband') 3907 3908 ! FIXME 3909 if (Dtset%symsigma/=0 .and. Sigp%nomegasr/=0) then 3910 if (cryst%idx_spatial_inversion() == 0) then 3911 write(msg,'(5a)')' setup_sigma : BUG :',ch10,& 3912 'It is not possible to use symsigma /= 0 to calculate the spectral function ',ch10,& 3913 'when the system does not have the spatial inversion. Please use symsigma=0 ' 3914 ABI_WARNING(msg) 3915 end if 3916 end if 3917 3918 if (mod10==SIG_GW_AC) then 3919 !if (Sigp%gwcalctyp/=1) ABI_ERROR("Self-consistency with AC not implemented") ! MRM: let's allow it 3920 if (Sigp%gwcomp==1) ABI_ERROR("AC with extrapolar technique not implemented") 3921 end if 3922 3923 call gaps%free() 3924 3925 ABI_FREE(val_indices) 3926 3927 DBG_EXIT('COLL') 3928 3929 end subroutine setup_sigma
ABINIT/setup_vcp [ Functions ]
NAME
setup_vcp
FUNCTION
Initialize the Vcp and compute the corresponding residue.
INPUTS
Dtset<type(dataset_type)>=all input variables for this dataset Gsph_c<gsphere_t>=Info on the G-sphere for W and Sigma_c Gsph_x<gsphere_t>=Info on the G-sphere for and Sigma_x Kmesh <kmesh_t>=Structure describing the k-point sampling. Qmesh <kmesh_t>=Structure describing the q-point sampling. Cryst<crystal_t>=Info on unit cell and symmetries. comm=Information about the xmpi_world
OUTPUT
Vcp_full<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique. Vcp_ks<vcoul_t>= Datatype gathering information on the coulombian interaction and the cutoff technique. coef_hyb=real variable containing the amount of GLOBAL hybridization
ABINIT/sigma_bksmask [ Functions ]
NAME
sigma_bksmask
FUNCTION
Compute tables for the distribution and the storage of the wavefunctions in the SIGMA code.
INPUTS
Dtset<type(dataset_type)>=all input variables for this dataset Sigp<sigparams_t>=Parameters governing the self-energy calculation. Kmesh <kmesh_t>=Structure describing the k-point sampling. nprocs=Total number of MPI processors my_rank=Rank of this this processor.
OUTPUT
my_spins(:)=Pointer to NULL in input. In output: list of spins treated by this node. bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will treat this state. keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)=True if this node will store this state in real space. ierr=Exit status.
SOURCE
4019 subroutine sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr) 4020 4021 !Arguments ------------------------------------ 4022 !scalars 4023 integer,intent(in) :: my_rank,nprocs 4024 integer,intent(out) :: ierr 4025 type(Dataset_type),intent(in) :: Dtset 4026 type(sigparams_t),intent(in) :: Sigp 4027 type(kmesh_t),intent(in) :: Kmesh 4028 !arrays 4029 integer,allocatable,intent(out) :: my_spins(:) 4030 logical,intent(out) :: bks_mask(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol) 4031 logical,intent(out) :: keep_ur(Sigp%nbnds,Kmesh%nibz,Sigp%nsppol) 4032 4033 !Local variables------------------------------- 4034 !scalars 4035 integer :: my_nspins,my_maxb,my_minb,isp,spin,band,nsppol,rank_spin 4036 logical :: store_ur 4037 !arrays 4038 integer :: tmp_spins(Sigp%nsppol),nprocs_spin(Sigp%nsppol) 4039 4040 ! ************************************************************************* 4041 4042 ierr=0; nsppol=Sigp%nsppol 4043 4044 ! List of spins for each node, number of processors per each spin 4045 ! and the MPI rank in the "spin" communicator. 4046 my_nspins=nsppol 4047 ABI_MALLOC(my_spins, (nsppol)) 4048 my_spins= [(isp, isp=1,nsppol)] 4049 nprocs_spin = nprocs; rank_spin = my_rank 4050 4051 if (nsppol==2 .and. nprocs>1) then 4052 ! Distribute spins (optimal distribution if nprocs is even) 4053 nprocs_spin(1) = nprocs/2 4054 nprocs_spin(2) = nprocs - nprocs/2 4055 my_nspins=1 4056 my_spins(1)=1 4057 if (my_rank+1>nprocs/2) then 4058 ! I will treate spin=2, compute shifted rank. 4059 my_spins(1)=2 4060 rank_spin = my_rank - nprocs/2 4061 end if 4062 end if 4063 4064 store_ur = (MODULO(Dtset%gwmem,10)==1) 4065 keep_ur=.FALSE.; bks_mask=.FALSE. 4066 4067 select case (Dtset%gwpara) 4068 case (1) 4069 ! Parallelization over transitions **without** memory distributions (Except for the spin). 4070 my_minb=1; my_maxb=Sigp%nbnds 4071 if (dtset%ucrpa>0) then 4072 bks_mask(my_minb:my_maxb,:,:)=.TRUE. 4073 if (store_ur) keep_ur(my_minb:my_maxb,:,:)=.TRUE. 4074 else 4075 do isp=1,my_nspins 4076 spin = my_spins(isp) 4077 bks_mask(my_minb:my_maxb,:,spin)=.TRUE. 4078 if (store_ur) keep_ur(my_minb:my_maxb,:,spin)=.TRUE. 4079 end do 4080 end if 4081 4082 case (2) 4083 ! Distribute bands and spin (alternating planes of bands) 4084 do isp=1,my_nspins 4085 spin = my_spins(isp) 4086 do band=1,Sigp%nbnds 4087 if (xmpi_distrib_with_replicas(band,Sigp%nbnds,rank_spin,nprocs_spin(spin))) then 4088 !if (MODULO(band-1,nprocs_spin(spin))==rank_spin) then 4089 bks_mask(band,:,spin)=.TRUE. 4090 if (store_ur) keep_ur(band,:,spin)=.TRUE. 4091 end if 4092 end do 4093 end do 4094 4095 #if 0 4096 ! Each node store the full set of occupied states to speed up Sigma_x. 4097 do isp=1,my_nspins 4098 spin = my_spins(isp) 4099 do ik_ibz=1,Kmesh%nibz 4100 ks_iv=ks_vbik(ik_ibz,spin) ! Valence index for this (k,s) 4101 bks_mask(1:ks_iv,:,spin)=.TRUE. 4102 if (store_ur) keep_ur(1:ks_iv,:,spin)=.TRUE. 4103 end do 4104 end do 4105 #endif 4106 4107 case default 4108 ABI_WARNING("Wrong value for gwpara") 4109 ierr = 1 4110 end select 4111 4112 ! Return my_spins with correct size. 4113 tmp_spins(1:my_nspins) = my_spins(1:my_nspins) 4114 4115 ABI_FREE(my_spins) 4116 ABI_MALLOC(my_spins, (my_nspins)) 4117 my_spins = tmp_spins(1:my_nspins) 4118 4119 end subroutine sigma_bksmask
ABINIT/sigma_tables [ Functions ]
NAME
sigma_tables
FUNCTION
Build symmetry tables used to speedup self-consistent GW calculations.
INPUTS
Kmesh <kmesh_t>=Structure describing the k-point sampling. [esymm(Kmesh%nibz,Sigp%nsppol)]= Bands_Symmetries SiDE EFFECTS Sigp<sigparams_t>=This routine initializes the tables: %Sigcij_tab %Sigxij_tab that are used to select the matrix elements of the self-energy that have to be calculated.
SOURCE
3953 subroutine sigma_tables(Sigp, Kmesh, esymm) 3954 3955 use m_sigtk, only : sigtk_sigma_tables 3956 3957 !Arguments ------------------------------------ 3958 !scalars 3959 type(sigparams_t),target,intent(inout) :: Sigp 3960 type(kmesh_t),intent(in) :: Kmesh 3961 !arrays 3962 type(esymm_t),optional,intent(in) :: esymm(Kmesh%nibz, Sigp%nsppol) 3963 3964 !Local variables------------------------------- 3965 !scalars 3966 integer :: ikcalc,ik_ibz 3967 logical :: sigc_is_herm, only_diago 3968 !arrays 3969 integer,allocatable :: kcalc2ibz(:) 3970 3971 ! ************************************************************************* 3972 3973 only_diago = sigp%gwcalctyp < 20 3974 sigc_is_herm = sigma_is_herm(Sigp) 3975 3976 ABI_MALLOC(kcalc2ibz, (sigp%nkptgw)) 3977 do ikcalc=1,sigp%nkptgw 3978 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 3979 kcalc2ibz(ikcalc) = ik_ibz 3980 end do 3981 3982 if (present(esymm)) then 3983 call sigtk_sigma_tables(sigp%nkptgw, kmesh%nibz, sigp%nsppol, sigp%minbnd, sigp%maxbnd, kcalc2ibz, & 3984 only_diago, sigc_is_herm, sigp%sigxij_tab, sigp%sigcij_tab, esymm=esymm) 3985 else 3986 call sigtk_sigma_tables(sigp%nkptgw, kmesh%nibz, sigp%nsppol, sigp%minbnd, sigp%maxbnd, kcalc2ibz, & 3987 only_diago, sigc_is_herm, sigp%sigxij_tab, sigp%sigcij_tab) 3988 end if 3989 3990 ABI_FREE(kcalc2ibz) 3991 3992 end subroutine sigma_tables
m_sigma_driver/sigma [ Functions ]
[ Top ] [ m_sigma_driver ] [ Functions ]
NAME
sigma
FUNCTION
Calculate the matrix elements of the self-energy operator.
INPUTS
acell(3)=length scales of primitive translations (bohr) codvsn=code version Dtfil<type(datafiles_type)>=variables related to files Dtset<type(dataset_type)>=all input variables for this dataset Pawang<type(pawang_type)>=paw angular mesh and related data Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data Psps<type(pseudopotential_type)>=variables related to pseudopotentials Before entering the first time in sigma, a significant part of Psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp. All the remaining components of Psps are to be initialized in the call to pspini. The next time the code enters screening, Psps might be identical to the one of the previous Dtset, in which case, no reinitialisation is scheduled in pspini.F90. rprim(3,3)=dimensionless real space primitive translations
OUTPUT
Output is written on the main abinit output file. Some results are stored in external files
NOTES
ON THE USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
171 subroutine sigma(acell,codvsn,Dtfil,Dtset,Pawang,Pawrad,Pawtab,Psps,rprim) 172 173 !Arguments ------------------------------------ 174 !scalars 175 character(len=8),intent(in) :: codvsn 176 type(Datafiles_type),intent(in) :: Dtfil 177 type(Dataset_type),intent(inout) :: Dtset 178 type(Pawang_type),intent(inout) :: Pawang 179 type(Pseudopotential_type),intent(inout) :: Psps 180 !arrays 181 real(dp),intent(in) :: acell(3),rprim(3,3) 182 type(Pawrad_type),intent(inout) :: Pawrad(Psps%ntypat*Psps%usepaw) 183 type(Pawtab_type),intent(inout) :: Pawtab(Psps%ntypat*Psps%usepaw) 184 185 !Local variables------------------------------- 186 !scalars 187 integer,parameter :: tim_fourdp5 = 5, master = 0, cplex1 = 1, ipert0 = 0, idir0 = 0, optrhoij1 = 1, ndat1 = 1 188 integer :: approx_type,b1gw,b2gw,cplex,cplex_dij,cplex_rhoij !,band 189 integer :: dim_kxcg,gwcalctyp,gnt_option,has_dijU,has_dijso,iab,bmin,bmax,irr_idx1,irr_idx2 190 integer :: iat,ib,ib1,ib2,ic,id_required,ider,ii,ik,ierr,ount 191 integer :: ik_bz,ikcalc,ik_ibz,ikxc,npw_k,omp_ncpus,pwx,ibz 192 integer :: isp,is_idx,istep,itypat,itypatcor,izero,jj,first_band,last_band 193 integer :: ks_iv,lcor,lmn2_size_max,mband,my_nband 194 integer :: mgfftf,mod10,moved_atm_inside,moved_rhor,n3xccc !,mgfft 195 integer :: nbsc,ndij,ndim,nfftf,nfftf_tot,nkcalc,gwc_nfft,gwc_nfftot,gwx_nfft,gwx_nfftot 196 integer :: ngrvdw,nhatgrdim,nkxc,nkxc1,nprocs,nscf,nspden_rhoij,nzlmopt,optene 197 integer :: optcut,optgr0,optgr1,optgr2,option,option_test,option_dij,optrad,psp_gencond 198 integer :: my_rank,rhoxsp_method,comm,use_aerhor,use_umklp,usexcnhat 199 integer :: ioe0j,spin,io,jb,nomega_sigc 200 integer :: temp_unt,ncid 201 integer :: work_size,nstates_per_proc,my_nbks 202 !integer :: jb_qp,ib_ks,ks_irr 203 integer :: gw1rdm,x1rdm 204 real(dp) :: compch_fft,compch_sph,r_s,rhoav,alpha 205 real(dp) :: drude_plsmf,my_plsmf,ecore,ecut_eff,ecutdg_eff,ehartree 206 real(dp) :: etot_sd,etot_mbb,evextnl_energy,ex_energy,gsqcutc_eff,gsqcutf_eff,gsqcut_shp,norm,oldefermi 207 real(dp) :: eh_energy,ekin_energy,evext_energy,den_int,coef_hyb,exc_mbb_energy,tol_empty 208 real(dp) :: ucvol,vxcavg,vxcavg_qp 209 real(dp) :: gwc_gsq,gwx_gsq,gw_gsq, gsqcut,boxcut,ecutf 210 real(dp) :: eff,mempercpu_mb,max_wfsmem_mb,nonscal_mem,ug_mem,ur_mem,cprj_mem 211 complex(dpc) :: max_degw,cdummy 212 logical :: rdm_update,readchkprdm,prtchkprdm 213 logical :: use_paw_aeur,dbg_mode,pole_screening,call_pawinit,is_dfpt=.false. 214 character(len=500) :: msg 215 character(len=fnlen) :: wfk_fname,pawden_fname,gw1rdm_fname 216 type(kmesh_t) :: Kmesh,Qmesh 217 type(ebands_t) :: ks_ebands, qp_ebands 218 type(vcoul_t) :: Vcp, Vcp_ks, Vcp_full 219 type(crystal_t) :: Cryst 220 type(Energies_type) :: KS_energies,QP_energies 221 type(Epsilonm1_results) :: Er 222 type(gsphere_t) :: Gsph_Max,Gsph_x,Gsph_c 223 type(hdr_type) :: Hdr_wfk,Hdr_sigma,Hdr_rhor 224 type(melflags_t) :: KS_mflags,QP_mflags 225 type(melements_t) :: KS_me, QP_me, GW1RDM_me 226 type(MPI_type) :: MPI_enreg_seq 227 type(paw_dmft_type) :: Paw_dmft 228 type(pawfgr_type) :: Pawfgr 229 type(ppmodel_t) :: PPm 230 type(sigparams_t) :: Sigp 231 type(sigma_t) :: Sr 232 type(wfdgw_t),target :: Wfd,Wfdf,Wfd_nato_master 233 type(wfdgw_t),pointer :: Wfd_nato_all 234 type(wave_t),pointer :: wave 235 type(wvl_data) :: Wvl 236 !arrays 237 integer :: gwc_ngfft(18),ngfftc(18),ngfftf(18),gwx_ngfft(18), units(2) 238 integer,allocatable :: sigmak_todo(:) 239 integer,allocatable :: nq_spl(:),nlmn_atm(:),my_spins(:) 240 integer,allocatable :: tmp_gfft(:,:),ks_vbik(:,:),nband(:,:),l_size_atm(:),qp_vbik(:,:) 241 integer,allocatable :: tmp_kstab(:,:,:),ks_irreptab(:,:,:),qp_irreptab(:,:,:),my_band_list(:) 242 real(dp),parameter :: k0(3) = zero 243 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2) 244 real(dp),allocatable :: weights(:),nat_occs(:,:),gw_rhor(:,:),gw_rhog(:,:),gw_vhartr(:) 245 real(dp),allocatable :: grchempottn(:,:),grewtn(:,:),grvdw(:,:),qmax(:) 246 real(dp),allocatable :: ks_nhat(:,:),ks_nhatgr(:,:,:),ks_rhog(:,:) 247 real(dp),allocatable :: ks_rhor(:,:),ks_vhartr(:),ks_vtrial(:,:),ks_vxc(:,:), ks_taur(:,:) 248 real(dp),allocatable :: kxc(:,:),qp_kxc(:,:),ph1d(:,:),ph1df(:,:) 249 real(dp),allocatable :: prev_rhor(:,:),prev_taur(:,:),qp_nhat(:,:) 250 real(dp),allocatable :: qp_nhatgr(:,:,:),qp_rhog(:,:),qp_rhor_paw(:,:) 251 real(dp),allocatable :: qp_rhor_n_one(:,:),qp_rhor_nt_one(:,:) 252 real(dp),allocatable :: qp_rhor(:,:),qp_vhartr(:),qp_vtrial(:,:),qp_vxc(:,:) 253 real(dp),allocatable :: qp_taur(:,:),igwene(:,:,:) 254 real(dp),allocatable :: vpsp(:),xccc3d(:),dijexc_core(:,:,:),dij_hf(:,:,:) 255 real(dp),allocatable :: nl_bks(:,:,:) 256 !real(dp),allocatable :: osoc_bks(:, :, :) 257 real(dp),allocatable :: ks_aepaw_rhor(:,:) !,ks_n_one_rhor(:,:),ks_nt_one_rhor(:,:) 258 complex(dpc) :: ovlp(2) 259 complex(dpc),allocatable :: ctmp(:,:),hbare(:,:,:,:) 260 complex(dpc),target,allocatable :: sigcme(:,:,:,:,:) 261 complex(dpc),allocatable :: hdft(:,:,:,:),htmp(:,:,:,:),uks2qp(:,:) 262 complex(dpc),allocatable :: xrdm_k_full(:,:,:), rdm_k(:,:), pot_k(:,:), nateigv(:,:,:,:), old_ks_purex(:,:), new_hartr(:,:) 263 complex(gwpc),allocatable :: kxcg(:,:),fxc_ADA(:,:,:) 264 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:) 265 complex(dpc),allocatable :: sigcme_k(:,:,:,:), rhot1_q_m(:,:,:,:,:,:,:), M1_q_m(:,:,:,:,:,:,:) 266 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:),bmask(:), bdm_mask(:,:,:),bdm2_mask(:,:,:) 267 type(esymm_t),target,allocatable :: KS_sym(:,:) 268 type(esymm_t),pointer :: QP_sym(:,:) 269 type(pawcprj_type),allocatable :: Cp1(:,:) !,Cp2(:,:) 270 type(littlegroup_t),allocatable :: Ltg_k(:) 271 type(Paw_an_type),allocatable :: KS_paw_an(:),QP_paw_an(:) 272 type(Paw_ij_type),allocatable :: KS_paw_ij(:),QP_paw_ij(:) 273 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 274 type(Pawrhoij_type),allocatable :: KS_Pawrhoij(:),QP_pawrhoij(:),prev_Pawrhoij(:),tmp_pawrhoij(:) 275 type(pawpwff_t),allocatable :: Paw_pwff(:) 276 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 277 type(plowannier_type) :: wanbz,wanibz,wanibz_in 278 type(operwan_realspace_type), allocatable :: rhot1(:,:) 279 280 !************************************************************************ 281 282 DBG_ENTER('COLL') 283 284 call timab(401,1,tsec) ! sigma(Total) 285 call timab(402,1,tsec) ! sigma(Init1) 286 units = [std_out, ab_out] 287 288 write(msg,'(7a)')& 289 ' SIGMA: Calculation of the GW corrections ',ch10,ch10,& 290 ' Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining.',ch10,& 291 ' Incorporated in ABINIT by V. Olevano, G.-M. Rignanese, and M. Torrent.' 292 call wrtout(units, msg) 293 294 if(dtset%ucrpa>0) then 295 write(msg,'(6a)')ch10,' cRPA Calculation: Calculation of the screened Coulomb interaction (ucrpa/=0) ',ch10 296 call wrtout(units, msg) 297 end if 298 299 #if defined HAVE_GW_DPC 300 if (gwpc /= 8) then 301 write(msg,'(6a)')ch10,& 302 'Number of bytes for double precision complex /=8 ',ch10,& 303 'Cannot continue due to kind mismatch in BLAS library ',ch10,& 304 'Some BLAS interfaces are not generated by abilint ' 305 ABI_ERROR(msg) 306 end if 307 write(msg,'(a,i2,a)')'.Using double precision arithmetic ; gwpc = ',gwpc,ch10 308 #else 309 write(msg,'(a,i2,a)')'.Using single precision arithmetic ; gwpc = ',gwpc,ch10 310 #endif 311 call wrtout(units, msg) 312 313 tol_empty=0.01 ! Initialize the tolerance used to decide if a band is empty (passed to m_sigx.F90) 314 gwcalctyp=Dtset%gwcalctyp 315 gw1rdm=Dtset%gw1rdm ! Input variable to decide if updates to the 1-RDM must be performed 316 x1rdm=Dtset%x1rdm ! Input variable to use pure exchange correction on the 1-RDM ( Sigma_x - Vxc ) 317 rdm_update=(gwcalctyp==21 .and. gw1rdm>0) ! Input variable to decide whether to update GW density matrix 318 readchkprdm=(Dtset%irdchkprdm==1) ! Input variable to decide if checkpoint files must be read 319 prtchkprdm=(Dtset%prtchkprdm==1) ! Input variable to decide if checkpoint files must be written 320 321 mod10 =MOD(Dtset%gwcalctyp,10) 322 323 ! Perform some additional checks for hybrid functional calculations 324 if (mod10 == 5) then 325 !if (Dtset%ixc_sigma<0 .and. .not.libxc_functionals_check()) then 326 ! ABI_ERROR('Hybrid functional calculations require the compilation with LIBXC library') 327 !end if 328 !XG 20171116 : I do not agree with this condition, as one might like to do a one-shot hybrid functional calculation 329 !on top of a LDA/GGA calculation ... give the power (and risks) to the user ! 330 !if(gwcalctyp<10) then 331 ! ABI_ERROR('gwcalctyp requires the update of energies and/or wavefunctions when performing hybrid XC calculations') 332 !end if 333 if (Dtset%usepaw == 1) then 334 ABI_ERROR('PAW version of hybrid functional calculations is not implemented') 335 end if 336 end if 337 338 !=== Initialize MPI variables, and parallelization level === 339 ! gwpara: 1--> parallelism over k-points, 2--> parallelism over bands. 340 ! In case of gwpara==1 memory is not parallelized. 341 ! If gwpara==2, bands are divided among processors but each proc has all the states where GW corrections are required. 342 comm = xmpi_world; my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 343 344 if (my_rank == master) then 345 wfk_fname = dtfil%fnamewffk 346 if (nctk_try_fort_or_ncfile(wfk_fname, msg) /= 0) then 347 ABI_ERROR(msg) 348 end if 349 end if 350 call xmpi_bcast(wfk_fname, master, comm, ierr) 351 352 ! === Some variables need to be initialized/nullify at start === 353 usexcnhat = 0 354 call energies_init(KS_energies) 355 call mkrdim(acell, rprim, rprimd) 356 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 357 ! 358 ! === Define FFT grid(s) sizes === 359 ! Be careful! This mesh is only used for densities, potentials and the matrix elements of v_Hxc. It is NOT the 360 ! (usually coarser) GW FFT mesh employed for the oscillator matrix elements that is defined in setmesh.F90. 361 ! See also NOTES in the comments at the beginning of this file. 362 ! NOTE: This mesh is defined in invars2m using ecutwfn, in GW Dtset%ecut is forced to be equal to Dtset%ecutwfn. 363 364 call pawfgr_init(Pawfgr, Dtset, mgfftf, nfftf, ecut_eff, ecutdg_eff, ngfftc, ngfftf, & 365 gsqcutc_eff=gsqcutc_eff, gsqcutf_eff=gsqcutf_eff, gmet=gmet, k0=k0) 366 367 ! Fake MPI_type for the sequential part. 368 call initmpi_seq(MPI_enreg_seq) 369 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all') 370 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 371 372 call print_ngfft(ngfftf,header='Dense FFT mesh used for densities and potentials') 373 nfftf_tot=PRODUCT(ngfftf(1:3)) 374 375 ! =========================================== 376 ! === Open and read pseudopotential files === 377 ! =========================================== 378 call pspini(dtset, dtfil, ecore, psp_gencond, gsqcutc_eff, gsqcutf_eff, pawrad, pawtab, psps, cryst%rprimd, comm_mpi=comm) 379 380 call timab(402,2,tsec) ! Init1 381 ! 382 ! =============================================== 383 ! ==== Initialize Sigp, Er and basic objects ==== 384 ! =============================================== 385 ! * Sigp is completetly initialized here. 386 ! * Er is only initialized with dimensions, (SCR|SUSC) file is read in mkdump_Er 387 call timab(403,1,tsec) ! setup_sigma 388 389 call setup_sigma(codvsn,wfk_fname,acell,rprim,Dtset,Dtfil,Psps,Pawtab,& 390 gwx_ngfft,gwc_ngfft,Hdr_wfk,Hdr_sigma,Cryst,Kmesh,Qmesh,ks_ebands,Gsph_Max,Gsph_x,Gsph_c,Vcp,Er,Sigp,comm) 391 392 call timab(403,2,tsec) ! setup_sigma 393 call timab(402,1,tsec) ! Init1 394 395 if (nprocs > Sigp%nbnds) then 396 write(msg,"(2(a,i0))")"The number of MPI procs: ", nprocs, " is greater than nband: ", sigp%nbnds 397 ABI_ERROR(msg) 398 end if 399 400 pole_screening = .FALSE. 401 if (Er%fform==2002) then 402 pole_screening = .TRUE. 403 ABI_WARNING(' EXPERIMENTAL - Using a pole-fit screening!') 404 end if 405 406 call print_ngfft(gwc_ngfft,header='FFT mesh for oscillator strengths used for Sigma_c') 407 call print_ngfft(gwx_ngfft,header='FFT mesh for oscillator strengths used for Sigma_x') 408 409 b1gw = Sigp%minbdgw; b2gw = Sigp%maxbdgw 410 411 gwc_nfftot=PRODUCT(gwc_ngfft(1:3)) 412 gwc_nfft =gwc_nfftot !no FFT // 413 414 gwx_nfftot=PRODUCT(gwx_ngfft(1:3)) 415 gwx_nfft =gwx_nfftot !no FFT // 416 417 ! TRYING TO RECREATE AN "ABINIT ENVIRONMENT" 418 KS_energies%e_corepsp=ecore/Cryst%ucvol 419 420 ! === Calculate KS occupation numbers and ks_vbk(nkibz,nsppol) ==== 421 ! * ks_vbk gives the (valence|last Fermi band) index for each k and spin. 422 ! * spinmagntarget is passed to fermi.F90 to fix the problem with newocc in case of magnetic metals 423 ABI_MALLOC(ks_vbik, (ks_ebands%nkpt, ks_ebands%nsppol)) 424 ABI_MALLOC(qp_vbik, (ks_ebands%nkpt, ks_ebands%nsppol)) 425 426 !call ebands_update_occ(ks_ebands,Dtset%spinmagntarget,prtvol=0) 427 ks_vbik(:,:) = ebands_get_valence_idx(ks_ebands) 428 429 ! ============================ 430 ! ==== PAW initialization ==== 431 ! ============================ 432 if (dtset%usepaw == 1) then 433 call chkpawovlp(cryst%natom, cryst%ntypat, dtset%pawovlp, pawtab, cryst%rmet, cryst%typat, cryst%xred) 434 435 cplex_dij = dtset%nspinor; cplex = 1; ndij = 1 436 437 ABI_MALLOC(ks_pawrhoij, (cryst%natom)) 438 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij, nspden_rhoij=nspden_rhoij, & 439 nspden=dtset%nspden, spnorb=dtset%pawspnorb, cpxocc=dtset%pawcpxocc) 440 call pawrhoij_alloc(ks_pawrhoij, cplex_rhoij, nspden_rhoij, dtset%nspinor, dtset%nsppol, cryst%typat, pawtab=pawtab) 441 442 ! Test if we have to call pawinit 443 gnt_option = 1; if (dtset%pawxcdev == 2 .or. (dtset%pawxcdev == 1 .and. dtset%positron /= 0)) gnt_option = 2 444 call paw_gencond(dtset, gnt_option, "test", call_pawinit) 445 446 if (psp_gencond == 1 .or. call_pawinit) then 447 call timab(553, 1, tsec) 448 gsqcut_shp = two * abs(dtset%diecut) * dtset%dilatmx**2 / pi**2 449 call pawinit(dtset%effmass_free, gnt_option, gsqcut_shp, zero, dtset%pawlcutd, dtset%pawlmix, & 450 psps%mpsang, dtset%pawnphi, cryst%nsym, dtset%pawntheta, pawang, pawrad, & 451 dtset%pawspnorb, pawtab, dtset%pawxcdev, dtset%ixc, dtset%usepotzero) 452 call timab(553,2,tsec) 453 454 ! Update internal values 455 call paw_gencond(dtset, gnt_option, "save", call_pawinit) 456 else 457 if (pawtab(1)%has_kij ==1) pawtab(1:cryst%ntypat)%has_kij = 2 458 if (pawtab(1)%has_nabla==1) pawtab(1:cryst%ntypat)%has_nabla = 2 459 end if 460 461 psps%n1xccc = maxval(pawtab(1:cryst%ntypat)%usetcore) 462 463 ! Initialize optional flags in Pawtab to zero 464 ! Cannot be done in Pawinit since the routine is called only if some parts are changed 465 pawtab(:)%has_nabla = 0 466 pawtab(:)%lamb_shielding = zero 467 468 call setsym_ylm(gprimd, pawang%l_max-1, cryst%nsym, dtset%pawprtvol, cryst%rprimd, cryst%symrec, pawang%zarot) 469 470 ! Initialize and compute data for DFT+U 471 Paw_dmft%use_dmft=Dtset%usedmft 472 call pawpuxinit(Dtset%dmatpuopt,Dtset%exchmix,Dtset%f4of2_sla,Dtset%f6of2_sla,& 473 is_dfpt,Dtset%jpawu,Dtset%lexexch,Dtset%lpawu,dtset%nspinor,Cryst%ntypat,dtset%optdcmagpawu,Pawang,Dtset%pawprtvol,& 474 Pawrad,Pawtab,Dtset%upawu,Dtset%usedmft,Dtset%useexexch,Dtset%usepawu,dtset%ucrpa) 475 476 if (my_rank == master) call pawtab_print(Pawtab) 477 478 ! Get Pawrhoij from the header of the WFK file. 479 call pawrhoij_copy(Hdr_wfk%pawrhoij, KS_Pawrhoij) 480 481 ! Evaluate form factor of radial part of phi.phj - tphi.tphj. 482 ! The q-grid must contain the FFT mesh used for sigma_c and the G-sphere for the exchange part. 483 ! We use the FFT mesh for sigma_c since COHSEX and the extrapolar method require oscillator 484 ! strengths on the FFT mesh. 485 ABI_MALLOC(tmp_gfft,(3, gwc_nfftot)) 486 call get_gftt(gwc_ngfft, k0, gmet, gwc_gsq, tmp_gfft) 487 ABI_FREE(tmp_gfft) 488 489 ! Set up q-grid, make qmax 20% larger than largest expected. 490 ABI_MALLOC(nq_spl, (Psps%ntypat)) 491 ABI_MALLOC(qmax, (Psps%ntypat)) 492 gwx_gsq = Dtset%ecutsigx / (two*pi**2) 493 gw_gsq = max(gwx_gsq, gwc_gsq) 494 qmax = sqrt(gw_gsq)*1.2d0 495 nq_spl = Psps%mqgrid_ff 496 ! write(std_out,*)"using nq_spl",nq_spl,"qmax=",qmax 497 498 rhoxsp_method = 1 ! Arnaud-Alouani (default in sigma) 499 !rhoxsp_method = 2 ! Shiskin-Kresse 500 if (dtset%pawoptosc /= 0) rhoxsp_method = dtset%pawoptosc 501 502 ABI_MALLOC(paw_pwff, (psps%ntypat)) 503 call pawpwff_init(paw_pwff, rhoxsp_method, nq_spl, qmax, gmet, pawrad, pawtab, psps) 504 505 ABI_FREE(nq_spl) 506 ABI_FREE(qmax) 507 508 ! Variables/arrays related to the fine FFT grid 509 ABI_CALLOC(ks_nhat, (nfftf, Dtset%nspden)) 510 511 ABI_MALLOC(pawfgrtab, (cryst%natom)) 512 call pawtab_get_lsize(pawtab, l_size_atm, cryst%natom, cryst%typat) 513 514 cplex = 1 515 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Dtset%nspden,Dtset%typat) 516 ABI_FREE(l_size_atm) 517 compch_fft=greatest_real 518 usexcnhat = MAXVAL(Pawtab(:)%usexcnhat) 519 ! * 0 if Vloc in atomic data is Vbare (Blochl's formulation) 520 ! * 1 if Vloc in atomic data is VH(tnzc) (Kresse's formulation) 521 call wrtout(std_out, sjoin(' using usexcnhat: ', itoa(usexcnhat))) 522 ! 523 ! Identify parts of the rectangular grid where the density has to be calculated === 524 optcut = 0; optgr0 = Dtset%pawstgylm; optgr1 = 0; optgr2 = 0; optrad = 1 - Dtset%pawstgylm 525 if (Dtset%pawcross==1) optrad=1 526 if (Dtset%xclevel==2.and.usexcnhat>0) optgr1=Dtset%pawstgylm 527 528 call nhatgrid(Cryst%atindx1, gmet, Cryst%natom, Cryst%natom, Cryst%nattyp, ngfftf, Cryst%ntypat,& 529 optcut, optgr0, optgr1, optgr2, optrad, Pawfgrtab, Pawtab, Cryst%rprimd, Cryst%typat, Cryst%ucvol, Cryst%xred) 530 531 call pawfgrtab_print(Pawfgrtab,Cryst%natom,unit=std_out,prtvol=Dtset%pawprtvol) 532 533 else 534 ABI_MALLOC(Paw_pwff, (0)) 535 ABI_MALLOC(Pawfgrtab, (0)) 536 end if ! End of PAW Initialization 537 538 ! Consistency check and additional stuff done only for GW with PAW. 539 ABI_MALLOC(Paw_onsite, (0)) 540 541 if (Dtset%usepaw == 1) then 542 if (Dtset%ecutwfn < Dtset%ecut) then 543 write(msg,"(3a)")& 544 " It is highly recommended to use ecutwfn = ecut for GW calculations with PAW since ",ch10,& 545 " an excessive truncation of the planewave basis set can lead to unphysical results." 546 ABI_WARNING(msg) 547 end if 548 549 ABI_CHECK(Dtset%useexexch == 0, "LEXX not yet implemented in GW") 550 ABI_CHECK(Paw_dmft%use_dmft == 0, "DMFT + GW not available") 551 552 ! Optionally read core orbitals from file and calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for the HF decoupling. 553 if (Sigp%use_sigxcore == 1) then 554 lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 555 ABI_MALLOC(dijexc_core,(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) 556 557 call paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,Dtset%prtvol,Psps%filpsp) 558 end if ! HF decoupling 559 560 if (Dtset%pawcross==1) then 561 ABI_SFREE(Paw_onsite) 562 ABI_MALLOC(Paw_onsite,(Cryst%natom)) 563 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat, & 564 Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab) 565 end if 566 end if 567 568 ! Allocate these arrays anyway, since they are passed to subroutines. 569 if (.not.allocated(ks_nhat)) then 570 ABI_MALLOC(ks_nhat, (nfftf, 0)) 571 end if 572 if (.not.allocated(dijexc_core)) then 573 ABI_MALLOC(dijexc_core, (1, 1, 0)) 574 end if 575 576 ! ================================================== 577 ! ==== Read KS band structure from the KSS file ==== 578 ! ================================================== 579 ! 580 ! Initialize Wfd, allocate wavefunctions and precalculate tables to do the FFT using the coarse gwc_ngfft. 581 mband=Sigp%nbnds 582 ABI_MALLOC(bks_mask, (mband, Kmesh%nibz, Sigp%nsppol)) 583 ABI_MALLOC(keep_ur , (mband, Kmesh%nibz, Sigp%nsppol)) 584 keep_ur=.FALSE.; bks_mask=.FALSE. 585 586 if (rdm_update) then 587 ABI_MALLOC(bdm_mask, (mband, Kmesh%nibz, Sigp%nsppol)) 588 bdm_mask=.FALSE. 589 if (my_rank==master) then 590 bdm_mask=.TRUE. 591 end if 592 ABI_MALLOC(bdm2_mask, (mband,Kmesh%nibz,Sigp%nsppol)) 593 bdm2_mask=.FALSE. 594 end if 595 ABI_MALLOC(nband,(Kmesh%nibz, Sigp%nsppol)) 596 nband=mband 597 598 ! autoparal section 599 if (dtset%max_ncpus/=0) then 600 ount =ab_out 601 ! Temporary table needed to estimate memory 602 ABI_MALLOC(nlmn_atm,(Cryst%natom)) 603 if (Dtset%usepaw==1) then 604 do iat=1,Cryst%natom 605 nlmn_atm(iat)=Pawtab(Cryst%typat(iat))%lmn_size 606 end do 607 end if 608 609 write(ount,'(a)')"--- !Autoparal" 610 write(ount,"(a)")'#Autoparal section for Sigma runs.' 611 write(ount,"(a)") "info:" 612 write(ount,"(a,i0)")" autoparal: ",dtset%autoparal 613 write(ount,"(a,i0)")" max_ncpus: ",dtset%max_ncpus 614 write(ount,"(a,i0)")" gwpara: ",dtset%gwpara 615 write(ount,"(a,i0)")" nkpt: ",dtset%nkpt 616 write(ount,"(a,i0)")" nsppol: ",dtset%nsppol 617 write(ount,"(a,i0)")" nspinor: ",dtset%nspinor 618 write(ount,"(a,i0)")" nbnds: ",Sigp%nbnds 619 620 work_size = mband * Kmesh%nibz * Sigp%nsppol 621 622 ! Non-scalable memory in Mb i.e. memory that is not distribute with MPI. 623 nonscal_mem = (two*gwpc*Er%npwe**2*Er%nomega*(Er%mqmem+1)*b2Mb) * 1.1_dp 624 625 ! List of configurations. 626 ! Assuming an OpenMP implementation with perfect speedup! 627 write(ount,"(a)")"configurations:" 628 do ii=1,dtset%max_ncpus 629 nstates_per_proc = 0 630 eff = HUGE(one) 631 max_wfsmem_mb = zero 632 633 do my_rank=0,ii-1 634 call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,ii,my_spins,bks_mask,keep_ur,ierr) 635 ABI_FREE(my_spins) 636 if (ierr /= 0) exit 637 my_nbks = COUNT(bks_mask) 638 nstates_per_proc = MAX(nstates_per_proc, my_nbks) 639 eff = MIN(eff, (one * work_size) / (ii * nstates_per_proc)) 640 641 ! Memory needed for Fourier components ug. 642 ug_mem = two*gwpc*Dtset%nspinor*Sigp%npwwfn*my_nbks*b2Mb 643 ! Memory needed for real space ur (use gwc_nfft, instead of gwx_nfft) 644 ur_mem = two*gwpc*Dtset%nspinor*gwc_nfft*COUNT(keep_ur)*b2Mb 645 ! Memory needed for PAW projections Cprj 646 cprj_mem = zero 647 if (Dtset%usepaw==1) cprj_mem = dp*Dtset%nspinor*SUM(nlmn_atm)*my_nbks*b2Mb 648 max_wfsmem_mb = MAX(max_wfsmem_mb, ug_mem + ur_mem + cprj_mem) 649 end do 650 if (ierr /= 0) cycle 651 652 ! Add the non-scalable part and increase by 10% to account for other datastructures. 653 mempercpu_mb = (max_wfsmem_mb + nonscal_mem) * 1.1_dp 654 do omp_ncpus=1,xomp_get_max_threads() 655 write(ount,"(a,i0)")" - tot_ncpus: ",ii * omp_ncpus 656 write(ount,"(a,i0)")" mpi_ncpus: ",ii 657 write(ount,"(a,i0)")" omp_ncpus: ",omp_ncpus 658 write(ount,"(a,f12.9)")" efficiency: ",eff 659 write(ount,"(a,f12.2)")" mem_per_cpu: ",mempercpu_mb 660 end do 661 end do 662 write(ount,'(a)')"..." 663 664 ABI_FREE(nlmn_atm) 665 ABI_ERROR_NODUMP("aborting now") 666 667 else 668 call sigma_bksmask(Dtset,Sigp,Kmesh,my_rank,nprocs,my_spins,bks_mask,keep_ur,ierr) 669 ABI_CHECK(ierr==0, "Error in sigma_bksmask") 670 end if 671 672 ! Each core stores the wavefunctions where GW corrections are required. 673 do isp=1,SIZE(my_spins) 674 spin = my_spins(isp) 675 do ikcalc=1,Sigp%nkptgw 676 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW 677 ii=Sigp%minbnd(ikcalc,spin); jj=Sigp%maxbnd(ikcalc,spin) 678 bks_mask(ii:jj,ik_ibz,spin) = .TRUE. 679 if (MODULO(Dtset%gwmem,10)==1) keep_ur(ii:jj,ik_ibz,spin)=.TRUE. 680 end do 681 end do 682 683 ABI_FREE(my_spins) 684 685 call wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,& 686 Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,& 687 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm,use_fnl_dir0der0=(Dtset%gw1rdm==2)) 688 689 ! MRM: also initialize the Wfd_nato_master for GW 1-RDM if required. 690 ! Warning, this should be replaced by copy but copy fails due to bands being allocated in different manners. 691 ! FIXME: Do it in the future! 692 if (rdm_update) then 693 bdm2_mask=bks_mask ! As bks_mask is going to be removed, save it in bdm2_mask to use it in Evext_nl 694 call wfd_init(Wfd_nato_master,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bdm_mask,& 695 Dtset%nspden,Dtset%nspinor,Dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,& 696 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,xmpi_comm_self)!comm) ! MPI_COMM_SELF 697 call Wfd_nato_master%read_wfk(wfk_fname,iomode_from_fname(wfk_fname)) 698 end if 699 700 if (Dtset%pawcross==1) then 701 call wfd_init(Wfdf,Cryst,Pawtab,Psps,keep_ur,mband,nband,Kmesh%nibz,Sigp%nsppol,bks_mask,& 702 Dtset%nspden,Dtset%nspinor,dtset%ecutwfn,Dtset%ecutsm,Dtset%dilatmx,Hdr_wfk%istwfk,Kmesh%ibz,gwc_ngfft,& 703 Dtset%nloalg,Dtset%prtvol,Dtset%pawprtvol,comm) 704 end if 705 706 ABI_FREE(bks_mask) 707 ABI_FREE(nband) 708 ABI_FREE(keep_ur) 709 710 call timab(402,2,tsec) ! sigma(Init1) 711 call timab(404,1,tsec) ! rdkss 712 713 call wfd%read_wfk(wfk_fname, iomode_from_fname(wfk_fname)) 714 715 if (Dtset%pawcross==1) then 716 call wfdgw_copy(Wfd, Wfdf) 717 call wfdf%change_ngfft(Cryst,Psps,ngfftf) 718 end if 719 720 ! This test has been disabled (too expensive!) 721 if (.False.) call wfd%test_ortho(Cryst,Pawtab,unit=ab_out,mode_paral="COLL") 722 723 call timab(404,2,tsec) ! rdkss 724 call timab(405,1,tsec) ! Init2 725 726 ! ============================================================== 727 ! ==== Find little group of the k-points for GW corrections ==== 728 ! ============================================================== 729 ! * The little group is used only if symsigma == 1 730 ! * If use_umklp == 1 then symmetries requiring an umklapp to preserve k_gw are included as well. 731 ! 732 ABI_MALLOC(Ltg_k, (Sigp%nkptgw)) 733 use_umklp = 1 734 do ikcalc=1,Sigp%nkptgw 735 if (Sigp%symsigma /= 0) then 736 call Ltg_k(ikcalc)%init(Sigp%kptgw(:,ikcalc), Qmesh%nbz, Qmesh%bz, Cryst, use_umklp, npwe=0) 737 end if 738 end do 739 740 ! Compute structure factor phases and large sphere cut-off 741 ABI_MALLOC(ph1d,(2, 3 * (2 * Dtset%mgfft + 1) * Cryst%natom)) 742 ABI_MALLOC(ph1df,(2, 3 * (2 * mgfftf + 1) * Cryst%natom)) 743 744 call getph(Cryst%atindx,Cryst%natom,ngfftc(1),ngfftc(2),ngfftc(3),ph1d,Cryst%xred) 745 746 if (Psps%usepaw == 1.and. Pawfgr%usefinegrid == 1) then 747 call getph(Cryst%atindx,Cryst%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,Cryst%xred) 748 else 749 ph1df(:,:)=ph1d(:,:) 750 end if 751 752 !=================================================================================== 753 !==== Classify the GW wavefunctions according to the irreducible representation ==== 754 !=================================================================================== 755 !* Warning still under development. 756 !* Only for SCGW. 757 !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 758 759 ABI_MALLOC(KS_sym,(Wfd%nkibz,Wfd%nsppol)) 760 761 if (Sigp%symsigma==1.and.gwcalctyp>=20) then 762 ! call check_zarot(Gsph_c%ng,Cryst,gwc_ngfft,Gsph_c%gvec,Psps,Pawang,Gsph_c%rottb,Gsph_c%rottbm1) 763 use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric 764 do spin=1,Wfd%nsppol 765 do ikcalc=1,Sigp%nkptgw 766 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 767 first_band = Sigp%minbnd(ikcalc,spin) 768 last_band = Sigp%maxbnd(ikcalc,spin) 769 ! call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,ks_ebands,Pawtab,Pawrad,Pawang,Psps,& 770 call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,ks_ebands,Pawtab,Pawrad,Pawang,Psps,& 771 & Dtset%tolsym,KS_sym(ik_ibz,spin)) 772 end do 773 end do 774 ! Recreate the Sig_ij tables taking advantage of the classification of the bands. 775 call sigma_tables(Sigp,Kmesh, esymm=KS_sym) 776 end if 777 778 call timab(405,2,tsec) ! Init2 779 call timab(406,1,tsec) ! make_vhxc 780 781 !=========================== 782 !=== COMPUTE THE DENSITY === 783 !=========================== 784 ! Evaluate the planewave part (complete charge in case of NC pseudos). 785 ABI_MALLOC(ks_rhor, (nfftf, dtset%nspden)) 786 ABI_MALLOC(ks_taur, (nfftf, dtset%nspden * dtset%usekden)) 787 788 call wfd%mkrho(cryst, psps, ks_ebands, ngfftf, nfftf, ks_rhor) 789 if ((rdm_update .and. Dtset%prtden /= 0) .and. Wfd%my_rank == master) then 790 ! Print initial (KS) density file as read (usefull to compare DEN files, cubes, etc.) 791 gw1rdm_fname = trim(dtfil%fnameabo_ks_den) ! and used on Sigma grids 792 call fftdatar_write("density",gw1rdm_fname,dtset%iomode,hdr_sigma,& 793 Cryst,ngfftf,cplex1,nfftf,dtset%nspden,ks_rhor,mpi_enreg_seq,ebands=ks_ebands) 794 end if 795 796 if (Dtset%usekden == 1) call wfd%mkrho(cryst, psps, ks_ebands, ngfftf, nfftf, ks_taur, optcalc=1) 797 798 !======================================== 799 !==== Additional computation for PAW ==== 800 !======================================== 801 nhatgrdim = 0 802 if (Dtset%usepaw==1) then 803 ! Calculate the compensation charge nhat. 804 if (Dtset%xclevel==2) nhatgrdim = usexcnhat * Dtset%pawnhatxc 805 cplex = 1; ider = 2 * nhatgrdim; izero = 0 806 if (nhatgrdim > 0) then 807 ABI_MALLOC(ks_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim)) 808 end if 809 if (nhatgrdim == 0) then 810 ABI_MALLOC(ks_nhatgr,(0,0,0)) 811 end if 812 813 call pawmknhat(compch_fft,cplex,ider,idir0,ipert0,izero,Cryst%gprimd,& 814 Cryst%natom,Cryst%natom,nfftf,ngfftf,nhatgrdim,Dtset%nspden,Cryst%ntypat,Pawang,& 815 Pawfgrtab,ks_nhatgr,ks_nhat,KS_Pawrhoij,KS_Pawrhoij,Pawtab,k0,Cryst%rprimd,& 816 Cryst%ucvol,dtset%usewvl,Cryst%xred) 817 818 ! === Evaluate onsite energies, potentials, densities === 819 ! * Initialize variables/arrays related to the PAW spheres. 820 ! * Initialize also lmselect (index of non-zero LM-moments of densities). 821 ABI_MALLOC(KS_paw_ij, (Cryst%natom)) 822 has_dijso = Dtset%pawspnorb; has_dijU = merge(0, 1, Dtset%usepawu == 0) 823 824 call paw_ij_nullify(KS_paw_ij) 825 call paw_ij_init(KS_paw_ij,cplex,Dtset%nspinor,Dtset%nsppol,& 826 Dtset%nspden,Dtset%pawspnorb,Cryst%natom,Cryst%ntypat,Cryst%typat,Pawtab,& 827 has_dij=1,has_dijhartree=1,has_dijhat=1,has_dijxc=1,has_dijxc_hat=1,has_dijxc_val=1,& 828 has_dijso=has_dijso,has_dijU=has_dijU,has_exexch_pot=1,has_pawu_occ=1) 829 830 nkxc1 = 0 831 ABI_MALLOC(KS_paw_an, (Cryst%natom)) 832 call paw_an_nullify(KS_paw_an) 833 call paw_an_init(KS_paw_an,Cryst%natom,Cryst%ntypat,nkxc1,0,Dtset%nspden,& 834 cplex,Dtset%pawxcdev,Cryst%typat,Pawang,Pawtab,has_vxc=1,has_vxcval=1) 835 836 ! Calculate onsite vxc with and without core charge. 837 nzlmopt=-1; option=0; compch_sph=greatest_real 838 call pawdenpot(compch_sph,KS_energies%e_paw,KS_energies%e_pawdc,ipert0,& 839 Dtset%ixc,Cryst%natom,Cryst%natom,Dtset%nspden,& 840 Cryst%ntypat,Dtset%nucdipmom,nzlmopt,option,KS_Paw_an,KS_Paw_an,KS_paw_ij,& 841 Pawang,Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,& 842 Pawtab,Dtset%pawxcdev,Dtset%spnorbscl,Dtset%xclevel,Dtset%xc_denpos,Dtset%xc_taupos,& 843 Cryst%ucvol,Psps%znuclpsp) 844 845 else 846 ABI_MALLOC(ks_nhatgr, (0, 0, 0)) 847 ABI_MALLOC(KS_paw_ij, (0)) 848 ABI_MALLOC(KS_paw_an, (0)) 849 end if ! PAW 850 851 call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,ks_rhor,Cryst%ucvol,& 852 Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 853 854 ! For PAW, add the compensation charge on the FFT mesh, then get rho(G). 855 if (Dtset%usepaw==1) ks_rhor = ks_rhor + ks_nhat 856 857 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_rhor,ucvol=ucvol) 858 if (Dtset%usekden==1) call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,ks_taur,optrhor=1,ucvol=ucvol) 859 860 ! FFT n(r) --> n(g) 861 ABI_MALLOC(ks_rhog, (2, nfftf)) 862 call fourdp(1, ks_rhog, ks_rhor(:,1),-1, MPI_enreg_seq, nfftf, 1, ngfftf, tim_fourdp5) 863 864 ! The following steps have been gathered in the setvtr routine: 865 ! - get Ewald energy and Ewald forces 866 ! - compute local ionic pseudopotential vpsp 867 ! - eventually compute 3D core electron density xccc3d 868 ! - eventually compute vxc and vhartr 869 ! - set up ks_vtrial 870 ! 871 !******************************************************************* 872 !**** NOTE THAT HERE Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 873 !******************************************************************* 874 875 ngrvdw = 0 876 ABI_MALLOC(grvdw, (3, ngrvdw)) 877 ABI_MALLOC(grchempottn, (3, Cryst%natom)) 878 ABI_MALLOC(grewtn, (3, Cryst%natom)) 879 nkxc = 0 880 if (Dtset%nspden == 1) nkxc = 2 881 if (Dtset%nspden >= 2) nkxc = 3 ! check GGA and spinor, quite a messy part!!! 882 ! In case of MGGA, fxc and kxc are not available and we dont need them in sigma (for now ...) 883 if (Dtset%ixc < 0 .and. libxc_functionals_ismgga()) nkxc = 0 884 if (nkxc /= 0) then 885 ABI_MALLOC(kxc, (nfftf, nkxc)) 886 end if 887 888 n3xccc = 0; if (Psps%n1xccc /= 0) n3xccc = nfftf 889 ABI_MALLOC(xccc3d, (n3xccc)) 890 ABI_MALLOC(ks_vhartr, (nfftf)) 891 ABI_MALLOC(ks_vtrial, (nfftf, Dtset%nspden)) 892 ABI_MALLOC(vpsp, (nfftf)) 893 ABI_MALLOC(ks_vxc, (nfftf, Dtset%nspden)) 894 895 optene = 4; moved_atm_inside = 0; moved_rhor = 0; istep = 1 896 897 call setvtr(Cryst%atindx1,Dtset,KS_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 898 istep,kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 899 Cryst%nattyp,nfftf,ngfftf,ngrvdw,ks_nhat,ks_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 900 optene,Pawang,Pawrad,KS_Pawrhoij,Pawtab,ph1df,Psps,ks_rhog,ks_rhor,Cryst%rmet,Cryst%rprimd,strsxc,& 901 Cryst%ucvol,usexcnhat,ks_vhartr,vpsp,ks_vtrial,ks_vxc,vxcavg,Wvl,xccc3d,Cryst%xred,taur=ks_taur) 902 903 !============================ 904 !==== Compute KS PAW Dij ==== 905 !============================ 906 if (Dtset%usepaw == 1) then 907 call timab(561,1,tsec) 908 909 ! Calculate the unsymmetrized Dij. 910 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,& 911 Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 912 Dtset%nspden,Cryst%ntypat,KS_paw_an,KS_paw_ij,Pawang,Pawfgrtab,& 913 Dtset%pawprtvol,Pawrad,KS_Pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 914 k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),ks_vtrial,ks_vxc,Cryst%xred,& 915 nucdipmom=Dtset%nucdipmom) 916 917 ! Symmetrize KS Dij 918 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,& 919 Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,KS_paw_ij,Pawang,& 920 Dtset%pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 921 922 ! Output the pseudopotential strengths Dij and the augmentation occupancies Rhoij. 923 call pawprt(Dtset,Cryst%natom,KS_paw_ij,KS_Pawrhoij,Pawtab) 924 call timab(561,2,tsec) 925 end if 926 927 call timab(406,2,tsec) ! make_vhxc 928 929 !=== Calculate Vxc(b1,b2,k,s)=<b1,k,s|v_{xc}|b2,k,s> for all the states included in GW === 930 ! * This part is parallelized within wfd%comm since each node has all GW wavefunctions. 931 ! * Note that vH matrix elements are calculated using the true uncutted interaction. 932 call timab(407,1,tsec) ! vHxc_me 933 934 call KS_mflags%reset() 935 if (rdm_update) then 936 KS_mflags%has_hbare=1 937 KS_mflags%has_kinetic=1 938 end if 939 KS_mflags%has_vhartree=1 940 KS_mflags%has_vxc =1 941 KS_mflags%has_vxcval =1 942 if (Dtset%usepawu /= 0 ) KS_mflags%has_vu = 1 943 if (Dtset%useexexch /= 0) KS_mflags%has_lexexch= 1 944 if (Sigp%use_sigxcore == 1) KS_mflags%has_sxcore = 1 945 if (gwcalctyp<10 ) KS_mflags%only_diago = 1 ! off-diagonal elements only for SC on wavefunctions. 946 947 if (.FALSE.) then ! quick and dirty hack to test HF contribution. 948 ABI_WARNING("testing on-site HF") 949 lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 950 ABI_MALLOC(dij_hf,(cplex_dij*lmn2_size_max,ndij,Cryst%natom)) 951 call paw_dijhf(ndij,cplex_dij,1,lmn2_size_max,Cryst%natom,Cryst%ntypat,Pawtab,Pawrad,Pawang,& 952 KS_Pawrhoij,dij_hf,Dtset%prtvol) 953 954 do iat=1,Cryst%natom 955 itypat = Cryst%typat(iat) 956 ii = Pawtab(itypat)%lmn2_size 957 KS_Paw_ij(iat)%dijxc(:,:) = dij_hf(1:cplex_dij*ii,:,iat) 958 KS_Paw_ij(iat)%dijxc_hat(:,:) = zero 959 end do 960 ABI_FREE(dij_hf) 961 962 !option_dij=3 963 !call symdij(Cryst%gprimd,Cryst%indsym,ipert0,& 964 !& Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 965 !& KS_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,& 966 !& Cryst%symafm,Cryst%symrec) 967 end if 968 969 ABI_MALLOC(tmp_kstab, (2, Wfd%nkibz, Wfd%nsppol)) 970 tmp_kstab=0 971 do spin=1,Sigp%nsppol 972 do ikcalc=1,Sigp%nkptgw ! No spin dependent! 973 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 974 tmp_kstab(1, ik_ibz, spin) = Sigp%minbnd(ikcalc, spin) 975 tmp_kstab(2, ik_ibz, spin) = Sigp%maxbnd(ikcalc, spin) 976 end do 977 end do 978 979 call calc_vhxc_me(Wfd, KS_mflags, KS_me, Cryst, Dtset, nfftf, ngfftf, & 980 ks_vtrial, ks_vhartr, ks_vxc, Psps, Pawtab, KS_paw_an, Pawang, Pawfgrtab, KS_paw_ij, dijexc_core, & 981 ks_rhor, usexcnhat, ks_nhat, ks_nhatgr, nhatgrdim, tmp_kstab, taur=ks_taur) 982 ABI_FREE(tmp_kstab) 983 984 !#ifdef DEV_HAVE_SCGW_SYM 985 !Set KS matrix elements connecting different irreps to zero. Do not touch unknown bands!. 986 if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then 987 bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 988 ABI_MALLOC(ks_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol)) 989 ks_irreptab=0 990 do spin=1,Sigp%nsppol 991 do ikcalc=1,Sigp%nkptgw 992 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 993 first_band = Sigp%minbnd(ikcalc,spin) 994 last_band = Sigp%maxbnd(ikcalc,spin) 995 if (.not.esymm_failed(KS_sym(ik_ibz,spin))) then 996 ks_irreptab(first_band:last_band,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(first_band:last_band) 997 !ks_irreptab(bmin:bmax,ik_ibz,spin) = KS_sym(ik_ibz,spin)%b2irrep(bmin:bmax) 998 end if 999 end do 1000 end do 1001 call KS_me%zero(ks_irreptab) 1002 ABI_FREE(ks_irreptab) 1003 end if 1004 !#endif 1005 1006 call KS_me%print(header="Matrix elements in the KS basis set", prtvol=Dtset%prtvol) 1007 ! 1008 ! If possible, calculate the EXX energy between the frozen core 1009 ! and the valence electrons using KS wavefunctions 1010 ! MG: Be careful here, since ex_energy is meaningful only if all occupied states are calculated. 1011 if( KS_mflags%has_sxcore ==1 ) then 1012 ! TODO 1013 !ex_energy = mels_get_exene_core(KS_me,kmesh,ks_ebands) 1014 ex_energy=zero 1015 do spin=1,Sigp%nsppol 1016 do ik=1,Kmesh%nibz 1017 do ib=b1gw,b2gw 1018 if (Sigp%nsig_ab==1) then 1019 ex_energy = ex_energy + half*ks_ebands%occ(ib,ik,spin)*Kmesh%wt(ik)*KS_me%sxcore(ib,ib,ik,spin) 1020 else 1021 ex_energy = ex_energy + half*ks_ebands%occ(ib,ik,spin)*Kmesh%wt(ik)*SUM(KS_me%sxcore(ib,ib,ik,:)) 1022 end if 1023 end do 1024 end do 1025 end do 1026 write(msg,'(a,2(es16.6,a))')' CORE Exchange energy with KS wavefunctions: ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV' 1027 call wrtout(std_out, msg) 1028 end if 1029 1030 call timab(407,2,tsec) ! vHxc_me 1031 call timab(408,1,tsec) ! hqp_init 1032 1033 ! Do not break this coding! 1034 ! When gwcalctyp>10, the order of the bands can be interexchanged after 1035 ! the diagonalization. Therefore, we have to correctly assign the matrix elements to the corresponding 1036 ! bands and we cannot skip the following even though it looks unuseful. 1037 if (gwcalctyp >= 10) call wrtout(std_out, ch10//' *************** KS Energies *******************') 1038 1039 !=== qp_ebands stores energies and occ. used for the calculation === 1040 ! * Initialize qp_ebands with KS values. 1041 ! * In case of SC update qp_ebands using the QPS file. 1042 call ebands_copy(ks_ebands, qp_ebands) 1043 1044 ABI_MALLOC(qp_rhor, (nfftf, Dtset%nspden)) 1045 ABI_MALLOC(qp_taur, (nfftf, Dtset%nspden * Dtset%usekden)) 1046 QP_sym => KS_sym 1047 1048 if (gwcalctyp<10) then 1049 ! one-shot GW, just do a copy of the KS density. 1050 qp_rhor=ks_rhor 1051 if(Dtset%usekden==1)qp_taur=ks_taur 1052 QP_sym => KS_sym 1053 else 1054 ! Self-consistent GW. 1055 ! Read the unitary matrix and the QP energies of the previous step from the QPS file. 1056 call energies_init(QP_energies) 1057 QP_energies%e_corepsp=ecore/Cryst%ucvol 1058 1059 ! m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}> 1060 ! Initialize the QP amplitudes with KS wavefunctions. 1061 ABI_MALLOC(Sr%m_ks_to_qp, (Sigp%nbnds, Sigp%nbnds, Kmesh%nibz, Sigp%nsppol)) 1062 Sr%m_ks_to_qp=czero 1063 do ib=1,Sigp%nbnds 1064 Sr%m_ks_to_qp(ib,ib,:,:)=cone 1065 end do 1066 1067 ! Now read m_ks_to_qp and update the energies in qp_ebands. 1068 ! TODO switch on the renormalization of n in sigma. 1069 ABI_MALLOC(prev_rhor, (nfftf, Dtset%nspden)) 1070 ABI_MALLOC(prev_taur, (nfftf, Dtset%nspden*Dtset%usekden)) 1071 ABI_MALLOC(prev_Pawrhoij, (Cryst%natom*Psps%usepaw)) 1072 1073 call rdqps(qp_ebands,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,1,nscf,& 1074 nfftf,ngfftf,Cryst%ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,Sr%m_ks_to_qp,prev_rhor,prev_Pawrhoij) 1075 1076 !Find the irreps associated to the QP amplitudes starting from the analogous table for the KS states. 1077 !bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1078 !allocate(qp_irreptab(bmin:bmax,Kmesh%nibz,Sigp%nsppol)) 1079 !qp_irreptab=0 1080 !!qp_irreptab=ks_irreptab 1081 1082 !do jb_qp=bmin,bmax 1083 !do ib_ks=bmin,bmax 1084 !if (ABS(Sr%m_ks_to_qp(ib_ks,jb_qp,ik_ibz,spin)) > tol12) then ! jb_qp has same the same character as ib_ks. 1085 !ks_irr = ks_irreptab(ib_ks,ib_ks,ik_ibz,spin) 1086 !qp_irreptab(jb_qp,jb_qp,ik_ibz,spin) = ks_irr 1087 !do ii=bmin,bmax 1088 !if (ks_irr == ks_irreptab(ii,ib_ks,ik_ibz,spin)) then 1089 !qp_irreptab(jb_qp,ii,ik_ibz,spin) = ks_irr 1090 !end if 1091 !end do 1092 !end if 1093 !end do 1094 !end do 1095 1096 if (nscf==0) prev_rhor=ks_rhor 1097 if (nscf==0 .and. Dtset%usekden==1) prev_taur=ks_taur 1098 1099 if (nscf>0.and.gwcalctyp>=20.and. wfd%my_rank == master) then 1100 ! Print the unitary transformation on std_out. 1101 call show_QP(qp_ebands,Sr%m_ks_to_qp,fromb=Sigp%minbdgw,tob=Sigp%maxbdgw,unit=std_out,tolmat=0.001_dp) 1102 end if 1103 1104 ! Compute QP wfg as linear combination of KS states === 1105 ! * Wfd%ug is modified inside calc_wf_qp 1106 ! * For PAW, update also the on-site projections. 1107 ! * WARNING the first dimension of MPI_enreg MUST be Kmesh%nibz 1108 ! TODO here we should use nbsc instead of nbnds 1109 1110 call wfd%rotate(Cryst, Sr%m_ks_to_qp) 1111 1112 ! Reinit the storage mode of Wfd as ug have been changed === 1113 ! Update also the wavefunctions for GW corrections on each processor 1114 call wfd%reset_ur_cprj() 1115 1116 ! This test has been disabled (too expensive!) 1117 if (.False.) call wfd%test_ortho(Cryst, Pawtab, unit=std_out) 1118 1119 ! Compute QP occupation numbers. 1120 call wrtout(std_out,'sigma: calculating QP occupation numbers:') 1121 call ebands_update_occ(qp_ebands, Dtset%spinmagntarget, prtvol=0) 1122 qp_vbik(:,:) = ebands_get_valence_idx(qp_ebands) 1123 1124 ! #ifdef DEV_HAVE_SCGW_SYM 1125 ! Calculate the irreducible representations of the new QP amplitdues. 1126 if (Sigp%symsigma==1.and.gwcalctyp>=20) then 1127 ABI_MALLOC(QP_sym,(Wfd%nkibz,Wfd%nsppol)) 1128 use_paw_aeur=.FALSE. ! should pass ngfftf but the dense mesh is not forced to be symmetric 1129 do spin=1,Wfd%nsppol 1130 do ikcalc=1,Sigp%nkptgw 1131 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1132 ! Quick fix for SCGW+symm TODO fix properly! 1133 first_band = Sigp%minbnd(ikcalc,spin) 1134 last_band = Sigp%maxbnd(ikcalc,spin) 1135 ! first_band = MINVAL(Sigp%minbnd(:,spin)) 1136 ! last_band = MAXVAL(Sigp%maxbnd(:,spin)) 1137 ! call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,Cryst,qp_ebands,Pawtab,Pawrad,Pawang,Psps,& 1138 call classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,Wfd%ngfft,Cryst,qp_ebands,Pawtab,Pawrad,Pawang,Psps,& 1139 & Dtset%tolsym,QP_sym(ik_ibz,spin)) 1140 end do 1141 end do 1142 1143 ! Recreate the Sig_ij tables taking advantage of the classification of the bands. 1144 call sigma_tables(Sigp, Kmesh, esymm=QP_sym) 1145 end if 1146 ! #endif 1147 1148 ! Compute QP density using the updated wfg. 1149 call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, qp_rhor) 1150 if (Dtset%usekden==1) call wfd%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, qp_taur, optcalc=1) 1151 1152 ! ======================================== 1153 ! ==== QP self-consistent GW with PAW ==== 1154 ! ======================================== 1155 if (Dtset%usepaw==1) then 1156 ABI_MALLOC(qp_nhat,(nfftf,Dtset%nspden)) 1157 nhatgrdim=0; if (Dtset%xclevel==2) nhatgrdim=usexcnhat 1158 ABI_MALLOC(qp_nhatgr,(nfftf,Dtset%nspden,3*nhatgrdim)) 1159 1160 ABI_MALLOC(QP_pawrhoij,(Cryst%natom)) 1161 ABI_MALLOC(QP_paw_ij,(Cryst%natom)) 1162 ABI_MALLOC(QP_paw_an,(Cryst%natom)) 1163 1164 ! Calculate new QP quantities: nhat, nhatgr, rho_ij, paw_ij, and paw_an. 1165 call paw_qpscgw(Wfd,nscf,nfftf,ngfftf,Dtset,Cryst,Kmesh,Psps,qp_ebands,& 1166 Pawang,Pawrad,Pawtab,Pawfgrtab,prev_Pawrhoij,& 1167 QP_pawrhoij,QP_paw_ij,QP_paw_an,QP_energies,qp_nhat,nhatgrdim,qp_nhatgr,compch_sph,compch_fft) 1168 else 1169 ABI_MALLOC(qp_nhat, (0, 0)) 1170 ABI_MALLOC(qp_nhatgr, (0, 0, 0)) 1171 ABI_MALLOC(QP_pawrhoij, (0)) 1172 ABI_MALLOC(QP_paw_ij, (0)) 1173 ABI_MALLOC(QP_paw_an, (0)) 1174 end if 1175 1176 ! here I should renormalize the density 1177 call test_charge(nfftf,ks_ebands%nelect,Dtset%nspden,qp_rhor,Cryst%ucvol,& 1178 Dtset%usepaw,usexcnhat,Pawfgr%usefinegrid,compch_sph,compch_fft,drude_plsmf) 1179 1180 if (Dtset%usepaw==1) qp_rhor(:,:)=qp_rhor(:,:)+qp_nhat(:,:) ! Add the "hat" term. 1181 1182 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_rhor,ucvol=ucvol) 1183 if(Dtset%usekden==1) then 1184 call prtrhomxmn(std_out,MPI_enreg_seq,nfftf,ngfftf,Dtset%nspden,1,qp_taur,optrhor=1,ucvol=ucvol) 1185 end if 1186 1187 ! Simple mixing of the PW density to damp oscillations in the Hartree potential. 1188 if (nscf>0 .and. (ABS(Dtset%rhoqpmix-one)>tol12) ) then 1189 write(msg,'(2a,f6.3)')ch10,' sigma: mixing QP densities using rhoqpmix= ',Dtset%rhoqpmix 1190 call wrtout(std_out, msg) 1191 qp_rhor = prev_rhor + Dtset%rhoqpmix*(qp_rhor-prev_rhor) 1192 if(Dtset%usekden==1) qp_taur = prev_taur + Dtset%rhoqpmix*(qp_taur-prev_taur) ! mix taur. 1193 end if 1194 1195 ABI_FREE(prev_rhor) 1196 ABI_FREE(prev_taur) 1197 if (Psps%usepaw==1.and.nscf>0) then 1198 call pawrhoij_free(prev_pawrhoij) 1199 end if 1200 ABI_FREE(prev_pawrhoij) 1201 1202 ABI_MALLOC(qp_rhog,(2,nfftf)) 1203 call fourdp(1,qp_rhog,qp_rhor(:,1),-1,MPI_enreg_seq,nfftf,1,ngfftf,tim_fourdp5) 1204 1205 ! =========================================== 1206 ! ==== Optional output of the QP density ==== 1207 ! =========================================== 1208 if (Dtset%prtden/=0 .and. wfd%my_rank == master) then 1209 call fftdatar_write("qp_rhor",dtfil%fnameabo_qp_den,dtset%iomode,hdr_sigma,& 1210 cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor,mpi_enreg_seq,ebands=qp_ebands) 1211 end if 1212 1213 ! =========================================== 1214 ! === Optional output of the full QP density 1215 ! =========================================== 1216 if (Wfd%usepaw==1.and.Dtset%prtden==2) then 1217 ABI_MALLOC(qp_rhor_paw ,(nfftf,Wfd%nspden)) 1218 ABI_MALLOC(qp_rhor_n_one ,(nfftf,Wfd%nspden)) 1219 ABI_MALLOC(qp_rhor_nt_one,(nfftf,Wfd%nspden)) 1220 1221 call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,qp_nhat,& 1222 Wfd%nspinor,Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,QP_pawrhoij,Pawtab,Dtset%prtvol,& 1223 qp_rhor,qp_rhor_paw,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 1224 1225 ABI_FREE(qp_rhor_n_one) 1226 ABI_FREE(qp_rhor_nt_one) 1227 1228 if (Dtset%prtvol > 9) then 1229 ! Print a normalisation check 1230 norm = SUM(qp_rhor_paw(:,1))*Cryst%ucvol/PRODUCT(Pawfgr%ngfft(1:3)) 1231 write(msg,'(a,F8.4)') ' QUASIPARTICLE DENSITY CALCULATED - NORM OF DENSITY: ',norm 1232 call wrtout(std_out, msg) 1233 end if 1234 1235 ! Write the density to file 1236 if (my_rank==master) then 1237 call fftdatar_write("qp_pawrhor",dtfil%fnameabo_qp_pawden,dtset%iomode,hdr_sigma,& 1238 cryst,ngfftf,cplex1,nfftf,dtset%nspden,qp_rhor_paw,mpi_enreg_seq,ebands=qp_ebands) 1239 end if 1240 ABI_FREE(qp_rhor_paw) 1241 end if 1242 1243 nkxc=0 1244 if (Dtset%nspden==1) nkxc=2 1245 if (Dtset%nspden>=2) nkxc=3 !check GGA and spinor that is messy !!! 1246 !In case of MGGA, fxc and kxc are not available and we dont need them 1247 !for the screening part (for now ...) 1248 if (Dtset%ixc<0.and.libxc_functionals_ismgga()) nkxc=0 1249 if (nkxc/=0) then 1250 ABI_MALLOC(qp_kxc,(nfftf,nkxc)) 1251 end if 1252 1253 ! **** NOTE THAT Vxc CONTAINS THE CORE-DENSITY CONTRIBUTION **** 1254 n3xccc=0; if (Psps%n1xccc/=0) n3xccc=nfftf 1255 ABI_MALLOC(qp_vhartr,(nfftf)) 1256 ABI_MALLOC(qp_vtrial,(nfftf,Dtset%nspden)) 1257 ABI_MALLOC(qp_vxc,(nfftf,Dtset%nspden)) 1258 1259 optene=4; moved_atm_inside=0; moved_rhor=0; istep=1 1260 1261 call setvtr(Cryst%atindx1,Dtset,QP_energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcutf_eff,& 1262 istep,qp_kxc,mgfftf,moved_atm_inside,moved_rhor,MPI_enreg_seq,& 1263 Cryst%nattyp,nfftf,ngfftf,ngrvdw,qp_nhat,qp_nhatgr,nhatgrdim,nkxc,Cryst%ntypat,Psps%n1xccc,n3xccc,& 1264 optene,Pawang,Pawrad,QP_pawrhoij,Pawtab,ph1df,Psps,qp_rhog,qp_rhor,Cryst%rmet,& 1265 Cryst%rprimd,strsxc,Cryst%ucvol,usexcnhat,qp_vhartr,vpsp,qp_vtrial,qp_vxc,vxcavg_qp,Wvl,& 1266 xccc3d,Cryst%xred,taur=qp_taur) 1267 1268 ABI_SFREE(qp_kxc) 1269 1270 if (Dtset%usepaw==1) then 1271 call timab(561,1,tsec) 1272 1273 ! Compute QP Dij 1274 call pawdij(cplex,Dtset%enunit,Cryst%gprimd,ipert0,& 1275 Cryst%natom,Cryst%natom,nfftf,ngfftf(1)*ngfftf(2)*ngfftf(3),& 1276 Dtset%nspden,Cryst%ntypat,QP_paw_an,QP_paw_ij,Pawang,Pawfgrtab,& 1277 Dtset%pawprtvol,Pawrad,QP_pawrhoij,Dtset%pawspnorb,Pawtab,Dtset%pawxcdev,& 1278 k0,Dtset%spnorbscl,Cryst%ucvol,dtset%cellcharge(1),qp_vtrial,qp_vxc,Cryst%xred,& 1279 nucdipmom=Dtset%nucdipmom) 1280 1281 ! Symmetrize total Dij 1282 option_dij=0 1283 1284 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,& 1285 Cryst%natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,& 1286 QP_paw_ij,Pawang,Dtset%pawprtvol,Pawtab,Cryst%rprimd,& 1287 Cryst%symafm,Cryst%symrec) 1288 1289 ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij. 1290 call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab) 1291 call timab(561,2,tsec) 1292 end if 1293 1294 ehartree=half*SUM(qp_rhor(:,1)*qp_vhartr(:))/DBLE(nfftf)*Cryst%ucvol 1295 1296 write(msg,'(a,80a)')ch10,('-',ii=1,80) 1297 call wrtout(ab_out,msg) 1298 write(msg,'(5a,f9.4,3a,es21.14,2a,es21.14)')ch10,& 1299 ' QP results after the unitary transformation in the KS subspace: ',ch10,ch10,& 1300 ' Number of electrons = ',qp_rhog(1,1)*Cryst%ucvol,ch10,ch10,& 1301 ' QP Band energy [Ha] = ',ebands_get_bandenergy(qp_ebands),ch10,& 1302 ' QP Hartree energy [Ha] = ',ehartree 1303 call wrtout(ab_out,msg) 1304 write(msg,'(a,80a)')ch10,('-',ii=1,80) 1305 call wrtout(ab_out,msg) 1306 1307 ! TODO Since plasmonpole model 2-3-4 depend on the Fourier components of the density 1308 ! in case of self-consistency we might calculate here the ppm coefficients using qp_rhor 1309 end if ! gwcalctyp>=10 1310 1311 ! KS hamiltonian: hdft(b1,b1,k,s)= <b1,k,s|H_s|b1,k,s> 1312 ABI_MALLOC(hdft, (b1gw:b2gw, b1gw:b2gw, Kmesh%nibz, Sigp%nsppol*Sigp%nsig_ab)) 1313 hdft = czero 1314 1315 if (Dtset%nspinor == 1) then 1316 do spin=1,Sigp%nsppol 1317 do ik=1,Kmesh%nibz 1318 do ib=b1gw,b2gw 1319 hdft(ib, ib, ik, spin) = ks_ebands%eig(ib, ik, spin) 1320 end do 1321 end do 1322 end do 1323 else 1324 ! Spinorial case 1325 ! * Note that here vxc contains the contribution of the core. 1326 ! * Scale ovlp if orthonormalization is not satisfied as npwwfn might be < npwvec. 1327 ! TODO add spin-orbit case and gwpara 2 1328 ABI_MALLOC(my_band_list, (wfd%mband)) 1329 ABI_MALLOC(bmask, (wfd%mband)) 1330 bmask = .False.; bmask(b1gw:b2gw) = .True. 1331 1332 if (Wfd%usepaw == 1) then 1333 ABI_MALLOC(Cp1,(Wfd%natom, Wfd%nspinor)) 1334 call pawcprj_alloc(Cp1, 0, Wfd%nlmn_atm) 1335 end if 1336 1337 do spin=1,Sigp%nsppol 1338 do ik_ibz=1,Kmesh%nibz 1339 ! Distribute bands in [b1gw, b2gw] range 1340 call wfd%distribute_bands(ik_ibz, spin, my_nband, my_band_list, bmask=bmask) 1341 if (my_nband == 0) cycle 1342 npw_k = Wfd%npwarr(ik_ibz) 1343 do ii=1,my_nband ! ib=b1gw,b2gw in sequential 1344 ib = my_band_list(ii) 1345 ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave, msg) == 0, msg) 1346 ug1 => wave%ug 1347 cdummy = xdotc(npw_k*Wfd%nspinor,ug1,1,ug1,1) 1348 ovlp(1) = REAL(cdummy); ovlp(2) = AIMAG(cdummy) 1349 1350 if (Psps%usepaw==1) then 1351 call wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 1352 ovlp = ovlp + paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab) 1353 end if 1354 ! write(std_out,*)ovlp(1),ovlp(2) 1355 norm=DBLE(ovlp(1)+ovlp(2)) 1356 ovlp(1)=DBLE(ovlp(1)/norm); ovlp(2)=DBLE(ovlp(2)/norm) ! ovlp(2)=cone-ovlp(1) 1357 hdft(ib,ib,ik_ibz,1) = ks_ebands%eig(ib,ik_ibz,1)*ovlp(1) - KS_me%vxc(ib,ib,ik_ibz,3) 1358 hdft(ib,ib,ik_ibz,2) = ks_ebands%eig(ib,ik_ibz,1)*ovlp(2) - KS_me%vxc(ib,ib,ik_ibz,4) 1359 hdft(ib,ib,ik_ibz,3) = KS_me%vxc(ib,ib,ik_ibz,3) 1360 hdft(ib,ib,ik_ibz,4) = KS_me%vxc(ib,ib,ik_ibz,4) 1361 end do 1362 end do 1363 end do 1364 1365 call xmpi_sum(hdft, wfd%comm, ierr) 1366 1367 ABI_FREE(my_band_list) 1368 ABI_FREE(bmask) 1369 if (Wfd%usepaw==1) then 1370 call pawcprj_free(Cp1) 1371 ABI_FREE(Cp1) 1372 end if 1373 end if 1374 1375 ! Initialize Sigma results === 1376 ! TODO it is better if we use ragged arrays indexed by the k-point 1377 call sigma_init(Sigp,Kmesh%nibz,Dtset%usepawu,Sr) 1378 1379 ! Setup bare Hamiltonian := T + v_{loc} + v_{nl} + v_H. 1380 ! 1381 ! * The representation depends wheter we are updating the wfs or not. 1382 ! * ks_vUme is zero unless we are using DFT+U as starting point, see calc_vHxc_braket 1383 ! * Note that vH matrix elements are calculated using the true uncutted interaction 1384 ! This should be changed if the cutoff is also used in the GS run. 1385 1386 if (gwcalctyp < 10) then 1387 ! For one-shot GW use the KS representation. 1388 Sr%hhartree = hdft - KS_me%vxcval 1389 1390 ! Additional stuff for PAW 1391 ! * DFT +U Hamiltonian 1392 ! * LEXX. 1393 ! * Core contribution estimated using Fock exchange. 1394 if (Dtset%usepaw==1) then 1395 if (Sigp%use_sigxcore == 1) Sr%hhartree = hdft - (KS_me%vxc - KS_me%sxcore) 1396 if (Dtset%usepawu /= 0) Sr%hhartree = Sr%hhartree - KS_me%vu 1397 if (Dtset%useexexch /= 0) then 1398 ABI_ERROR("useexexch > 0 not implemented") 1399 Sr%hhartree = Sr%hhartree - KS_me%vlexx 1400 end if 1401 end if 1402 1403 else 1404 ! Self-consistent on energies and|or wavefunctions. 1405 ! * For NC get the bare Hamiltonian $H_{bare}= T+v_{loc}+ v_{nl}$ in the KS representation 1406 ! * For PAW, calculate the matrix elements of h0, store also the new Dij in QP_Paw_ij. 1407 ! * h0 is defined as T+vH[tn+nhat+tnZc] + vxc[tnc] + dij_eff and 1408 ! dij_eff = dij^0 + dij^hartree + dij^xc-dij^xc_val + dijhat - dijhat_val. 1409 ! In the above expression tn, tnhat are QP quantities. 1410 if (Dtset%usepaw==0) then 1411 ABI_MALLOC(hbare, (b1gw:b2gw,b1gw:b2gw,Kmesh%nibz,Sigp%nsppol*Sigp%nsig_ab)) 1412 hbare = hdft - KS_me%vhartree - KS_me%vxcval 1413 1414 ! Change basis from KS to QP, hbare is overwritten: A_{QP} = U^\dagger A_{KS} U 1415 ABI_MALLOC(htmp, (b1gw:b2gw,b1gw:b2gw, Kmesh%nibz, Sigp%nsppol*Sigp%nsig_ab)) 1416 ABI_MALLOC(ctmp, (b1gw:b2gw, b1gw:b2gw)) 1417 ABI_MALLOC(uks2qp, (b1gw:b2gw, b1gw:b2gw)) 1418 1419 htmp = hbare; hbare = czero 1420 1421 do spin=1,Sigp%nsppol 1422 do ik=1,Kmesh%nibz 1423 uks2qp(:,:) = Sr%m_ks_to_qp(b1gw:b2gw,b1gw:b2gw,ik,spin) 1424 do iab=1,Sigp%nsig_ab 1425 is_idx=spin; if (Sigp%nsig_ab>1) is_idx=iab 1426 ctmp(:,:)=MATMUL(htmp(:,:,ik,is_idx),uks2qp(:,:)) 1427 hbare(:,:,ik,is_idx)=MATMUL(TRANSPOSE(CONJG(uks2qp)),ctmp) 1428 end do 1429 end do 1430 end do 1431 ABI_FREE(htmp) 1432 ABI_FREE(ctmp) 1433 ABI_FREE(uks2qp) 1434 end if ! usepaw==0 1435 1436 ! Calculate the QP matrix elements 1437 ! This part is parallelized within MPI_COMM_WORD since each node has all GW wavefunctions. 1438 ! For PAW, we have to construct the new bare Hamiltonian. 1439 call wrtout(std_out,ch10//' *************** QP Energies *******************') 1440 1441 call QP_mflags%reset() 1442 ! if (gwcalctyp<20) QP_mflags%only_diago=1 ! For e-only, no need of off-diagonal elements. 1443 QP_mflags%has_vhartree=1 1444 if (Dtset%usepaw==1) then 1445 QP_mflags%has_kinetic =1 1446 QP_mflags%has_hbare =1 1447 end if 1448 !QP_mflags%has_vxc =1 1449 !QP_mflags%has_vxcval =1 1450 !if (Sigp%gwcalctyp >100) QP_mflags%has_vxcval_hybrid=1 1451 if (mod10==5 .and. & 1452 (Dtset%ixc_sigma==-402 .or. Dtset%ixc_sigma==-406 .or. Dtset%ixc_sigma==-427 .or. Dtset%ixc_sigma==-428 .or. & 1453 Dtset%ixc_sigma==-456 .or. Dtset%ixc_sigma==41 .or. Dtset%ixc_sigma==42)) then 1454 QP_mflags%has_vxcval_hybrid = 1 1455 end if 1456 1457 !if (Sigp%use_sigxcore==1) QP_mflags%has_sxcore =1 1458 !if (Dtset%usepawu/=0) QP_mflags%has_vu =1 1459 !if (Dtset%useexexch/=0) QP_mflags%has_lexexch=1 1460 1461 ABI_MALLOC(tmp_kstab, (2, Wfd%nkibz, Wfd%nsppol)) 1462 tmp_kstab=0 1463 do spin=1,Sigp%nsppol 1464 do ikcalc=1,Sigp%nkptgw ! No spin dependent! 1465 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1466 tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin) 1467 tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin) 1468 end do 1469 end do 1470 1471 call calc_vhxc_me(Wfd,QP_mflags,QP_me,Cryst,Dtset,nfftf,ngfftf,& 1472 qp_vtrial,qp_vhartr,qp_vxc,Psps,Pawtab,QP_paw_an,Pawang,Pawfgrtab,QP_paw_ij,dijexc_core,& 1473 qp_rhor,usexcnhat,qp_nhat,qp_nhatgr,nhatgrdim,tmp_kstab,taur=qp_taur) 1474 ABI_FREE(tmp_kstab) 1475 1476 ! #ifdef DEV_HAVE_SCGW_SYM 1477 if (gwcalctyp>=20 .and. Sigp%symsigma>0) then 1478 bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1479 ABI_MALLOC(qp_irreptab,(bmin:bmax,Kmesh%nibz,Sigp%nsppol)) 1480 qp_irreptab=0 1481 do spin=1,Sigp%nsppol 1482 do ikcalc=1,Sigp%nkptgw 1483 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1484 first_band = Sigp%minbnd(ikcalc,spin) 1485 last_band = Sigp%maxbnd(ikcalc,spin) 1486 if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then 1487 qp_irreptab(first_band:last_band,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(first_band:last_band) 1488 ! qp_irreptab(bmin:bmax,ik_ibz,spin) = QP_sym(ik_ibz,spin)%b2irrep(bmin:bmax) 1489 end if 1490 end do 1491 end do 1492 call QP_me%zero(qp_irreptab) 1493 ABI_FREE(qp_irreptab) 1494 end if 1495 ! #endif 1496 1497 call QP_me%print(header="Matrix elements in the QP basis set", prtvol=Dtset%prtvol) 1498 1499 ! Output the QP pseudopotential strengths Dij and the augmentation occupancies Rhoij. 1500 if (Dtset%usepaw==1) then 1501 call wrtout(std_out," *** After calc_vHxc_braket *** ") 1502 ! TODO terminate the implementation of this routine. 1503 call paw_ij_print(QP_Paw_ij,unit=std_out,pawprtvol=Dtset%pawprtvol,pawspnorb=Dtset%pawspnorb,mode_paral="COLL") 1504 call pawprt(Dtset,Cryst%natom,QP_paw_ij,QP_Pawrhoij,Pawtab) 1505 end if 1506 1507 if (Dtset%usepaw==0) then 1508 ! GA: We have an odd bug here. I have to unroll this loop, otherwise it 1509 ! might cause segfault when running on several nodes. 1510 ! 1511 ! Sr%hhartree = hbare + QP_me%vhartree 1512 if(QP_mflags%has_vxcval_hybrid==0) then 1513 do spin=1,Sigp%nsppol*Sr%nsig_ab 1514 do ikcalc=1,Sr%nkibz 1515 do ib1=b1gw,b2gw 1516 do ib2=b1gw,b2gw 1517 Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + QP_me%vhartree(ib2,ib1,ikcalc,spin) 1518 end do 1519 end do 1520 end do 1521 end do 1522 else 1523 do spin=1,Sigp%nsppol*Sr%nsig_ab 1524 do ikcalc=1,Sr%nkibz 1525 do ib1=b1gw,b2gw 1526 do ib2=b1gw,b2gw 1527 Sr%hhartree(ib2,ib1,ikcalc,spin) = hbare(ib2,ib1,ikcalc,spin) + & 1528 QP_me%vhartree(ib2,ib1,ikcalc,spin) + & 1529 QP_me%vxcval_hybrid(ib2,ib1,ikcalc,spin) 1530 end do 1531 end do 1532 end do 1533 end do 1534 end if 1535 else 1536 Sr%hhartree=QP_me%hbare 1537 end if 1538 1539 ! #ifdef DEV_HAVE_SCGW_SYM 1540 if (gwcalctyp>=20 .and. Sigp%symsigma > 0) then 1541 ! bmin=Sigp%minbdgw; bmax=Sigp%maxbdgw 1542 do spin=1,Sigp%nsppol 1543 do ik_ibz=1,Kmesh%nibz 1544 if (.not.esymm_failed(QP_sym(ik_ibz,spin))) then 1545 bmin=Sigp%minbnd(ik_ibz,spin); bmax=Sigp%minbnd(ik_ibz,spin) 1546 do ib2=bmin,bmax 1547 irr_idx2 = QP_sym(ik_ibz,spin)%b2irrep(ib2) 1548 do ib1=bmin,bmax 1549 irr_idx1 = QP_sym(ik_ibz,spin)%b2irrep(ib1) 1550 if (irr_idx1/=irr_idx2 .and. ALL((/irr_idx1,irr_idx2/)/=0) ) Sr%hhartree(ib1,ib2,ik_ibz,spin) = czero 1551 end do 1552 end do 1553 end if 1554 end do 1555 end do 1556 end if 1557 ! #endif 1558 1559 ABI_FREE(qp_rhog) 1560 ABI_FREE(qp_vhartr) 1561 ABI_FREE(qp_vxc) 1562 ABI_FREE(qp_nhat) 1563 ABI_FREE(qp_nhatgr) 1564 call QP_me%free() 1565 end if ! gwcalctyp<10 1566 1567 ! Free some memory 1568 ABI_SFREE(hbare) 1569 ABI_SFREE(hdft) 1570 1571 ! Prepare the storage of QP amplitudes and energies 1572 ! Initialize with KS wavefunctions and energies. 1573 Sr%eigvec_qp=czero; Sr%en_qp_diago=zero 1574 do ib=1,Sigp%nbnds 1575 Sr%en_qp_diago(ib,:,:) = ks_ebands%eig(ib,:,:) 1576 Sr%eigvec_qp(ib,ib,:,:) = cone 1577 end do 1578 1579 ! Store <n,k,s|V_xc[n_val]|n,k,s> and <n,k,s|V_U|n,k,s> === 1580 ! Note that we store the matrix elements of V_xc in the KS basis set, not in the QP basis set 1581 ! Matrix elements of V_U are zero unless we are using DFT+U as starting point 1582 do ib=b1gw,b2gw 1583 Sr%vxcme(ib,:,:)=KS_me%vxcval(ib,ib,:,:) 1584 if (Dtset%usepawu/=0) Sr%vUme (ib,:,:)=KS_me%vu(ib,ib,:,:) 1585 end do 1586 1587 ! Initial guess for the GW energies 1588 ! Save the energies of the previous iteration. 1589 do spin=1,Sigp%nsppol 1590 do ik=1,Kmesh%nibz 1591 do ib=1,Sigp%nbnds 1592 Sr%e0 (ib,ik,spin) = qp_ebands%eig(ib,ik,spin) 1593 Sr%egw(ib,ik,spin) = qp_ebands%eig(ib,ik,spin) 1594 end do 1595 Sr%e0gap(ik,spin) = zero 1596 ks_iv = ks_vbik(ik, spin) 1597 if (Sigp%nbnds>=ks_iv+1) Sr%e0gap(ik,spin)=Sr%e0(ks_iv+1,ik,spin)-Sr%e0(ks_iv,ik,spin) 1598 end do 1599 end do 1600 ! 1601 !=== If required apply a scissor operator or update the energies === 1602 !TODO check if other Sr entries have to be updated 1603 !moreover this part should be done only in case of semiconductors 1604 !FIXME To me it makes more sense if we apply the scissor to ks_ebands but I have to RECHECK csigme 1605 if (ABS(Sigp%mbpt_sciss)>tol6) then 1606 write(msg,'(6a,f10.5,2a)')ch10,& 1607 ' sigma : performing a first self-consistency',ch10,& 1608 ' update of the energies in G by a scissor operator',ch10, & 1609 ' applying a scissor operator of ',Sigp%mbpt_sciss*Ha_eV,' [eV] ',ch10 1610 call wrtout(units, msg) 1611 do spin=1,Sigp%nsppol 1612 do ik=1,Kmesh%nibz 1613 ks_iv=ks_vbik(ik,spin) 1614 if (Sigp%nbnds>=ks_iv+1) then 1615 Sr%egw (ks_iv+1:Sigp%nbnds,ik,spin) = Sr%egw (ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss 1616 qp_ebands%eig(ks_iv+1:Sigp%nbnds,ik,spin) = qp_ebands%eig(ks_iv+1:Sigp%nbnds,ik,spin)+Sigp%mbpt_sciss 1617 end if 1618 end do 1619 end do 1620 !call apply_scissor(qp_ebands,Sigp%mbpt_sciss) 1621 else if (.FALSE.) then 1622 write(msg,'(4a)')ch10,& 1623 ' sigma : performing a first self-consistency',ch10,& 1624 ' update of the energies in G by a previous GW calculation' 1625 call wrtout(units, msg) 1626 ! TODO Recheck this part, is not clear to me! 1627 ABI_MALLOC(igwene,(qp_ebands%mband, qp_ebands%nkpt, qp_ebands%nsppol)) 1628 call rdgw(qp_ebands, '__in.gw__', igwene, extrapolate=.TRUE.) 1629 ABI_FREE(igwene) 1630 Sr%egw=qp_ebands%eig 1631 ! 1632 ! * Recalculate the new fermi level. 1633 call ebands_update_occ(qp_ebands, Dtset%spinmagntarget, prtvol=0) 1634 end if 1635 1636 ! In case of AC refer all the energies wrt to the fermi level 1637 ! Be careful because results from ppmodel cannot be used for AC 1638 ! FIXME check ks_energy or qp_energy (in case of SCGW?) 1639 1640 if (mod10 == SIG_GW_AC) then 1641 ! All these quantities will be passed to csigme 1642 ! if I skipped the self-consistent part then here I have to use fermi 1643 qp_ebands%eig = qp_ebands%eig - qp_ebands%fermie 1644 Sr%egw = Sr%egw - qp_ebands%fermie 1645 Sr%e0 = Sr%e0 - qp_ebands%fermie 1646 oldefermi = qp_ebands%fermie 1647 ! TODO Recheck fermi 1648 ! Clean EVERYTHING in particulare the treatment of E fermi 1649 qp_ebands%fermie=zero 1650 end if 1651 ! 1652 ! === Setup frequencies around the KS\QP eigenvalues to compute Sigma derivatives (notice the spin) === 1653 ! TODO it is better using an odd Sr%nomega4sd so that the KS\QP eigenvalue is in the middle 1654 ioe0j=Sr%nomega4sd/2+1 1655 do spin=1,Sigp%nsppol 1656 do ikcalc=1,Sigp%nkptgw 1657 ib1=Sigp%minbnd(ikcalc,spin) 1658 ib2=Sigp%maxbnd(ikcalc,spin) 1659 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 1660 do jb=ib1,ib2 1661 do io=1,Sr%nomega4sd 1662 Sr%omega4sd(jb,ik_ibz,io,spin)=Sr%egw(jb,ik_ibz,spin)+Sigp%deltae*(io-ioe0j) 1663 end do 1664 end do 1665 end do 1666 end do 1667 1668 call timab(408,2,tsec) ! hqp_init 1669 call timab(409,1,tsec) ! getW 1670 1671 ! Get epsilon^{-1} either from the _SCR or the _SUSC file and store it in Er%epsm1 1672 ! If Er%mqmem==0, allocate and read a single q-slice inside csigme. 1673 ! TODO Er%nomega should be initialized so that only the frequencies really needed are stored in memory 1674 ! TODO The same piece of code is present in screening. 1675 if (sigma_needs_w(Sigp)) then 1676 select case (dtset%gwgamma) 1677 case (0) 1678 id_required=4; ikxc=0; approx_type=0; option_test=0; dim_kxcg=0 1679 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1680 1681 case (1, 2) 1682 ! ALDA TDDFT kernel vertex 1683 !ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1684 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1685 1686 if (Dtset%usepaw==1) then 1687 ! If we have PAW, we need the full density on the fine grid 1688 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1689 if (Dtset%getpawden==0 .and. Dtset%irdpawden==0) then 1690 ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1691 end if 1692 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin)) 1693 if (.not. file_exists(dtfil%filpawdensin)) then 1694 ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1695 end if 1696 1697 ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1698 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1699 ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm) 1700 1701 call Hdr_rhor%free() 1702 call pawrhoij_free(tmp_pawrhoij) 1703 ABI_FREE(tmp_pawrhoij) 1704 end if ! Dtset%usepaw==1 1705 1706 id_required=4; ikxc=7; approx_type=1; dim_kxcg=1 1707 if (Dtset%gwgamma==1) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1708 if (Dtset%gwgamma==2) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1709 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1710 1711 dbg_mode=.FALSE. 1712 if (Dtset%usepaw==1) then 1713 ! Use PAW all-electron density 1714 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_aepaw_rhor,& 1715 Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode) 1716 else 1717 ! Norm-conserving 1718 call kxc_driver(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,ks_rhor,& 1719 Er%npwe,dim_kxcg,kxcg,Gsph_c%gvec,xmpi_comm_self,dbg_mode=dbg_mode) 1720 end if 1721 1722 case (3, 4) 1723 ! ADA non-local kernel vertex 1724 !ABI_CHECK(Dtset%usepaw==0,"GWGamma + PAW not available") 1725 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1726 ABI_CHECK(Sigp%nsppol==1,"ADA vertex for GWGamma not available yet for spin-polarised cases") 1727 1728 if (Dtset%usepaw==1) then ! If we have PAW, we need the full density on the fine grid 1729 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Sigp%nsppol)) 1730 if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then 1731 ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1732 end if 1733 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin)) 1734 if (.not. file_exists(dtfil%filpawdensin)) then 1735 ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1736 end if 1737 1738 ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1739 1740 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1741 ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm) 1742 1743 call Hdr_rhor%free() 1744 call pawrhoij_free(tmp_pawrhoij) 1745 ABI_FREE(tmp_pawrhoij) 1746 end if ! Dtset%usepaw==1 1747 1748 id_required=4; ikxc=7; approx_type=2; 1749 if (Dtset%gwgamma==3) option_test=1 ! TESTELECTRON, vertex in chi0 *and* sigma 1750 if (Dtset%gwgamma==4) option_test=0 ! TESTPARTICLE, vertex in chi0 only 1751 ABI_MALLOC(fxc_ADA,(Er%npwe,Er%npwe,Er%nqibz)) 1752 ! Use userrd to set kappa 1753 if (Dtset%userrd==zero) Dtset%userrd = 2.1_dp 1754 ! Set correct value of kappa (should be scaled with alpha*r_s where) 1755 ! r_s is Wigner-Seitz radius and alpha=(4/(9*Pi))^(1/3) 1756 rhoav = (drude_plsmf*drude_plsmf)/four_pi 1757 r_s = (three/(four_pi*rhoav))**third 1758 alpha = (four*ninth*piinv)**third 1759 Dtset%userrd = Dtset%userrd/(alpha*r_s) 1760 1761 dbg_mode=.TRUE. 1762 if (Dtset%usepaw==1) then 1763 ! Use PAW all-electron density 1764 call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,& 1765 ks_aepaw_rhor,Er%npwe,Er%nqibz,Er%qibz,& 1766 fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode) 1767 else 1768 ! Norm conserving 1769 call kxc_ADA(Dtset,Cryst,ikxc,ngfftf,nfftf_tot,Wfd%nspden,& 1770 ks_rhor,Er%npwe,Er%nqibz,Er%qibz,& 1771 fxc_ADA,Gsph_c%gvec,xmpi_comm_self,kappa_init=Dtset%userrd,dbg_mode=dbg_mode) 1772 end if 1773 1774 dim_kxcg = 0 1775 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1776 1777 case (-3, -4, -5, -6, -7, -8) 1778 ! bootstrap kernels 1779 ! ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1780 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1781 1782 if (Dtset%usepaw==1) then 1783 ! If we have PAW, we need the full density on the fine grid 1784 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1785 if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then 1786 ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1787 end if 1788 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin)) 1789 if (.not. file_exists(dtfil%filpawdensin)) then 1790 ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1791 end if 1792 1793 ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1794 1795 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1796 ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm) 1797 1798 call Hdr_rhor%free() 1799 call pawrhoij_free(tmp_pawrhoij) 1800 ABI_FREE(tmp_pawrhoij) 1801 end if ! Dtset%usepaw==1 1802 1803 id_required=4; ikxc=7; dim_kxcg=0 1804 1805 if (dtset%gwgamma>-5) then 1806 approx_type=4 ! full fxc(G,G') 1807 else if (dtset%gwgamma>-7) then 1808 approx_type=5 ! fxc(0,0) one-shot 1809 else 1810 approx_type=6 ! rpa-type bootstrap 1811 end if 1812 1813 option_test=MOD(Dtset%gwgamma,2) 1814 ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma 1815 ! 0 -> TESTPARTICLE, vertex in chi0 only 1816 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1817 1818 case (-11) 1819 !LR+ALDA kernel 1820 !ABI_CHECK(Dtset%usepaw==0,"GWGamma=1 or 2 + PAW not available") 1821 ABI_CHECK(Er%ID==0,"Er%ID should be 0") 1822 1823 if (Dtset%usepaw==1) then 1824 ! If we have PAW, we need the full density on the fine grid 1825 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1826 if (Dtset%getpawden==0.and.Dtset%irdpawden==0) then 1827 ABI_ERROR("Must use get/irdpawden to provide a _PAWDEN file!") 1828 end if 1829 call wrtout(std_out,sjoin('Checking for existence of: ',Dtfil%filpawdensin)) 1830 if (.not. file_exists(dtfil%filpawdensin)) then 1831 ABI_ERROR(sjoin("Missing file:", dtfil%filpawdensin)) 1832 end if 1833 1834 ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1835 1836 call read_rhor(Dtfil%filpawdensin, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1837 ks_aepaw_rhor, hdr_rhor, tmp_pawrhoij, wfd%comm) 1838 1839 call hdr_rhor%free() 1840 call pawrhoij_free(tmp_pawrhoij) 1841 ABI_FREE(tmp_pawrhoij) 1842 end if ! Dtset%usepaw==1 1843 1844 id_required=4; ikxc=7; dim_kxcg=1 1845 approx_type=7 1846 option_test=1 1847 ! 1 -> TESTELECTRON, vertex in chi0 *and* sigma 1848 ABI_MALLOC(kxcg,(nfftf_tot,dim_kxcg)) 1849 1850 case default 1851 ABI_ERROR(sjoin("Wrong gwgamma:", itoa(dtset%gwgamma))) 1852 end select 1853 1854 ! Set plasma frequency 1855 my_plsmf = drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf = Dtset%ppmfrq 1856 Dtset%ppmfrq = my_plsmf 1857 1858 ! if(dtset%ucrpa==0) then 1859 if (Dtset%gwgamma<3) then 1860 call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,& 1861 approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,& 1862 nfftf_tot,ngfftf,comm) 1863 else 1864 call mkdump_Er(Er,Vcp,Er%npwe,Gsph_c%gvec,dim_kxcg,kxcg,id_required,& 1865 approx_type,ikxc,option_test,Dtfil%fnameabo_scr,Dtset%iomode,& 1866 nfftf_tot,ngfftf,comm,fxc_ADA) 1867 end if 1868 ! end if 1869 ABI_SFREE(kxcg) 1870 ABI_SFREE(fxc_ADA) 1871 ABI_SFREE(ks_aepaw_rhor) 1872 end if 1873 ! 1874 ! ================================================ 1875 ! ==== Calculate plasmonpole model parameters ==== 1876 ! ================================================ 1877 ! TODO Maybe its better if we use mqmem as input variable 1878 use_aerhor=0 1879 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden*use_aerhor)) 1880 1881 if (sigma_needs_ppm(Sigp)) then 1882 my_plsmf=drude_plsmf; if (Dtset%ppmfrq>tol6) my_plsmf=Dtset%ppmfrq 1883 call ppm_init(PPm,Er%mqmem,Er%nqibz,Er%npwe,Sigp%ppmodel,my_plsmf,Dtset%gw_invalid_freq) 1884 1885 ! PPm%force_plsmf= force_ppmfrq ! this line to change the plasme frequency in HL expression. 1886 1887 if (Wfd%usepaw==1 .and. Ppm%userho==1) then 1888 ! * For PAW and ppmodel 2-3-4 we need the AE rho(G) without compensation charge. 1889 ! * It would be possible to calculate rho(G) using Paw_pwff, though. It should be faster but 1890 ! results will depend on the expression used for the matrix elements. This approach is safer. 1891 use_aerhor=1 1892 ABI_FREE(ks_aepaw_rhor) 1893 ABI_MALLOC(ks_aepaw_rhor,(nfftf,Wfd%nspden)) 1894 1895 ! Check if input density file is available, otherwise compute 1896 pawden_fname = strcat(Dtfil%filnam_ds(3), '_PAWDEN') 1897 call wrtout(std_out,sjoin('Checking for existence of:',pawden_fname)) 1898 if (file_exists(pawden_fname)) then 1899 ! Read density from file 1900 ABI_MALLOC(tmp_pawrhoij,(cryst%natom*wfd%usepaw)) 1901 1902 call read_rhor(pawden_fname, cplex1, nfftf_tot, Wfd%nspden, ngfftf, 1, MPI_enreg_seq, & 1903 ks_aepaw_rhor, Hdr_rhor, tmp_pawrhoij, wfd%comm) 1904 1905 call Hdr_rhor%free() 1906 call pawrhoij_free(tmp_pawrhoij) 1907 ABI_FREE(tmp_pawrhoij) 1908 else 1909 ! Have to calculate PAW AW rhor from scratch 1910 ABI_MALLOC(qp_rhor_n_one,(pawfgr%nfft,Dtset%nspden)) 1911 ABI_MALLOC(qp_rhor_nt_one,(pawfgr%nfft,Dtset%nspden)) 1912 1913 ! FIXME 1914 ABI_WARNING(" denfgr in sigma seems to produce wrong results") 1915 write(std_out,*)" input tilde ks_rhor integrates: ",SUM(ks_rhor(:,1))*Cryst%ucvol/nfftf 1916 1917 call denfgr(Cryst%atindx1,Cryst%gmet,comm,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,ks_nhat,Wfd%nspinor,& 1918 Wfd%nsppol,Wfd%nspden,Cryst%ntypat,Pawfgr,Pawrad,KS_Pawrhoij,Pawtab,Dtset%prtvol, & 1919 ks_rhor,ks_aepaw_rhor,qp_rhor_n_one,qp_rhor_nt_one,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 1920 1921 ABI_FREE(qp_rhor_n_one) 1922 ABI_FREE(qp_rhor_nt_one) 1923 end if 1924 1925 write(msg,'(a,f8.4)')' sigma: PAW AE density used for PPmodel integrates to: ',SUM(ks_aepaw_rhor(:,1))*Cryst%ucvol/nfftf 1926 call wrtout(std_out, msg) 1927 1928 if (Er%mqmem/=0) then 1929 ! Calculate ppmodel parameters for all q-points. 1930 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_aepaw_rhor(:,1)) 1931 end if 1932 1933 else 1934 ! NC or PAW with PPmodel 1. 1935 if (Er%mqmem/=0) then 1936 ! Calculate ppmodel parameters for all q-points 1937 call setup_ppmodel(PPm,Cryst,Qmesh,Er%npwe,Er%nomega,Er%omega,Er%epsm1,nfftf,Gsph_c%gvec,ngfftf,ks_rhor(:,1)) 1938 end if 1939 end if ! PAW or NC PPm and/or needs density 1940 end if ! sigma_needs_ppm 1941 1942 call timab(409,2,tsec) ! getW 1943 1944 if (wfd%my_rank == master) then 1945 ! Write info on the run on ab_out, then open files to store final results. 1946 call ebands_report_gap(ks_ebands,header='KS Band Gaps',unit=ab_out) 1947 if(dtset%ucrpa == 0) then 1948 call write_sigma_header(Sigp,Er,Cryst,Kmesh,Qmesh) 1949 end if 1950 1951 ! unt_gw: File with GW corrections. 1952 ! unt_sig: Self-energy as a function of frequency. 1953 ! unt_sgr: Derivative wrt omega of the Self-energy. 1954 ! unt_sigc: Sigma_c(eik) MRM 1955 ! unt_sgm: Sigma on the Matsubara axis (imag axis) 1956 1957 if (open_file(Dtfil%fnameabo_gw, msg, unit=unt_gw, status='unknown', form='formatted') /= 0) then 1958 ABI_ERROR(msg) 1959 end if 1960 write(unt_gw,"(a)")"# QP energies E in eV" 1961 write(unt_gw,"(a)")"# Format:" 1962 write(unt_gw,"(a)")"# kpoint" 1963 write(unt_gw,"(a)")"# number of bands computed" 1964 write(unt_gw,"(2a)")"# band index, Re(E), E-E0, Im(E)", ch10 1965 write(unt_gw,"(2(i0,1x),a)")Sigp%nkptgw,Sigp%nsppol, "# nkptgw, nsppol" 1966 1967 if (open_file(Dtfil%fnameabo_gwdiag, msg, unit=unt_gwdiag, status='unknown', form='formatted') /= 0) then 1968 ABI_ERROR(msg) 1969 end if 1970 write(unt_gwdiag,*)Sigp%nkptgw,Sigp%nsppol 1971 1972 if (open_file(Dtfil%fnameabo_sig, msg, unit=unt_sig, status='unknown', form='formatted') /= 0) then 1973 ABI_ERROR(msg) 1974 end if 1975 write(unt_sig,"(a)")"# Sigma_xc and spectral function A along the real frequency axis in eV units" 1976 write(unt_sig,"(a)")"# Format:" 1977 write(unt_sig,"(a)")"# kpoint" 1978 write(unt_sig,"(a)")"# min_band(k) max_band(k)" 1979 write(unt_sig,"(a)")"# For each frequency w:" 1980 write(unt_sig,"(2a)")"# w, {Re(Sigma_b(w)), Im(Sigma_b(w), A_b(w) for b in [min_band, max_band]}",ch10 1981 1982 if (open_file(Dtfil%fnameabo_sgr, msg, unit=unt_sgr, status='unknown', form='formatted') /= 0) then 1983 ABI_ERROR(msg) 1984 end if 1985 write(unt_sgr,"(a)")"# Derivatives of Sigma_c(omega) wrt omega in eV units" 1986 1987 ! Sigma_c(eik) MRM 1988 if (open_file(trim(Dtfil%fnameabo_sgr)//'_SIGC', msg, unit=unt_sigc, status='unknown', form='formatted') /= 0) then 1989 ABI_ERROR(msg) 1990 end if 1991 1992 if (mod10 == SIG_GW_AC) then 1993 ! Sigma along the imaginary axis. 1994 if (open_file(Dtfil%fnameabo_sgm, msg, unit=unt_sgm, status='unknown', form='formatted') /= 0) then 1995 ABI_ERROR(msg) 1996 end if 1997 write(unt_sgm,"(a)")"# Sigma_xc along the imaginary frequency axis in eV units" 1998 end if 1999 end if 2000 2001 !======================================================================= 2002 !==== Calculate self-energy and output the results for each k-point ==== 2003 !======================================================================= 2004 ! Here it would be possible to calculate the QP correction for the same k-point using a PPmodel 2005 ! in the first iteration just to improve the initial guess and CD or AC in the second step. Useful 2006 ! if the KS level is far from the expected QP results. Previously it was possible to do such 2007 ! calculation by simply listing the same k-point twice in kptgw. Now this trick is not allowed anymore. 2008 ! Everything, indeed, should be done in a clean and transparent way inside csigme. 2009 2010 call wfd%print(mode_paral='PERS') 2011 2012 call wrtout(std_out,sigma_type_from_key(mod10)) 2013 2014 if (gwcalctyp<10) then 2015 msg = " Perturbative Calculation" 2016 if (gwcalctyp==1) msg = " Newton Raphson method " 2017 else if (gwcalctyp<20) then 2018 msg = " Self-Consistent on Energies only" 2019 else if (gwcalctyp==21) then 2020 msg = " GW 1-RDM Corrections" 2021 else 2022 msg = " Self-Consistent on Energies and Wavefunctions" 2023 end if 2024 if (Dtset%ucrpa==0) call wrtout(std_out,msg) 2025 2026 !================================================= 2027 !==== Calculate the matrix elements of Sigma ===== 2028 !================================================= 2029 2030 nomega_sigc=Sr%nomega_r+Sr%nomega4sd; if (mod10==SIG_GW_AC) nomega_sigc=Sr%nomega_i 2031 2032 ! min and max band indices for GW corrections. 2033 ib1=Sigp%minbdgw; ib2=Sigp%maxbdgw 2034 2035 !MG TODO: I don't like the fact that ib1 and ib2 are redefined here because this 2036 ! prevents me from refactoring the code. In particular I want to store the self-energy 2037 ! results inside the sigma_results datatypes hence one needs to know all the dimensions 2038 ! at the beginning of the execution (e.g. in setup_sigma) so that one can easily allocate the arrays in the type. 2039 if(Dtset%ucrpa>=1) then 2040 !Read the band 2041 if (dtset%plowan_compute<10)then 2042 if (open_file("forlb.ovlp",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then 2043 ABI_ERROR(msg) 2044 end if 2045 rewind(temp_unt) 2046 read(temp_unt,*) 2047 read(temp_unt,*) 2048 read(temp_unt,*) msg, ib1, ib2 2049 close(temp_unt) 2050 else 2051 if (open_file("data.plowann",msg,newunit=temp_unt,form="formatted", status="unknown") /= 0) then 2052 ABI_ERROR(msg) 2053 end if 2054 rewind(temp_unt) 2055 read(temp_unt,*) 2056 read(temp_unt,*) 2057 read(temp_unt,'(a7,2i4)') msg, ib1, ib2 2058 close(temp_unt) 2059 endif 2060 endif 2061 2062 ! Do not store it for rdm_update because it might be too large! 2063 if(.not.rdm_update) then 2064 ABI_MALLOC(sigcme,(nomega_sigc,ib1:ib2,ib1:ib2,Sigp%nkptgw,Sigp%nsppol*Sigp%nsig_ab)) 2065 sigcme=czero 2066 endif 2067 2068 ! if (.False. .and. psps%usepaw == 0 .and. wfd%nspinor == 1 .and. any(dtset%so_psp /= 0)) then 2069 ! call wrtout(std_out, "Computing SOC contribution with first-order perturbation theory") 2070 ! ABI_MALLOC(bks_mask, (wfd%mband, wfd%nkibz, wfd%nsppol)) 2071 ! bks_mask = .False. 2072 ! do spin=1,wfd%nsppol 2073 ! do ikcalc=1,Sigp%nkptgw 2074 ! ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW 2075 ! ii=Sigp%minbnd(ikcalc, spin); jj=Sigp%maxbnd(ikcalc, spin) 2076 ! bks_mask(ii:jj, ik_ibz, spin) = .True. 2077 ! end do 2078 ! end do 2079 ! 2080 ! call wfd_get_socpert(wfd, cryst, psps, pawtab, bks_mask, osoc_bks) 2081 ! ABI_FREE(bks_mask) 2082 ! ABI_FREE(osoc_bks) 2083 ! end if 2084 2085 !========================================================== 2086 !==== Exchange part using the dense gwx_ngfft FFT mesh ==== 2087 !========================================================== 2088 !TODO : load distribution is not optimal if band parallelism is used. 2089 !but this problem was also affecting the old implementation. 2090 call timab(421,1,tsec) ! calc_sigx_me 2091 2092 ! This routine shows how the wavefunctions are distributed. 2093 !call wfd_show_bkstab(Wfd,unit=std_out) 2094 2095 if(Dtset%ucrpa>=1) then 2096 ! Calculation of the Rho_n for calc_ucrpa 2097 call wrtout(std_out, sjoin("begin of Ucrpa calc for a nb of kpoints of: ",itoa(Sigp%nkptgw))) 2098 ! Wannier basis: rhot1_q_m will need an atom index in the near future. 2099 ic=0 2100 do itypat=1,cryst%ntypat 2101 if(pawtab(itypat)%lpawu.ne.-1) then 2102 ndim=2*dtset%lpawu(itypat)+1 2103 ic=ic+1 2104 itypatcor=itypat 2105 lcor=dtset%lpawu(itypat) 2106 end if 2107 end do 2108 if(ic>1) then 2109 ABI_ERROR("number of correlated species is larger than one") 2110 end if 2111 ABI_MALLOC(rhot1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz)) 2112 ABI_MALLOC(M1_q_m,(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,ndim,ndim,sigp%npwx,Qmesh%nibz)) 2113 2114 M1_q_m=czero 2115 rhot1_q_m=czero 2116 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 2117 !Initialization of wan objects and getting psichies 2118 !Allocation of rhot1 2119 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 2120 if (dtset%plowan_compute>=10)then 2121 call wrtout(units, " cRPA calculations using wannier weights from data.plowann") 2122 call init_plowannier(dtset%plowan_bandf,dtset%plowan_bandi,dtset%plowan_compute,dtset%plowan_iatom,& 2123 dtset%plowan_it,dtset%plowan_lcalc,dtset%plowan_natom,dtset%plowan_nbl,dtset%plowan_nt,& 2124 dtset%plowan_projcalc,dtset%acell_orig,dtset%kptns,dtset%nimage,dtset%nkpt,dtset%nspinor,& 2125 dtset%nsppol,dtset%wtk,dtset%dmft_t2g,wanibz_in) 2126 call get_plowannier(wanibz_in,wanibz,dtset) 2127 call fullbz_plowannier(dtset,kmesh,cryst,pawang,wanibz,wanbz) 2128 ABI_MALLOC(rhot1,(sigp%npwx,Qmesh%nibz)) 2129 do pwx=1,sigp%npwx 2130 do ibz=1,Qmesh%nibz 2131 call init_operwan_realspace(wanbz,rhot1(pwx,ibz)) 2132 call zero_operwan_realspace(wanbz,rhot1(pwx,ibz)) 2133 end do 2134 end do 2135 endif 2136 2137 ! do ikcalc=1,Sigp%nkptgw 2138 2139 ! if(cryst%nsym==1) nkcalc=Kmesh%nbz 2140 ! if(cryst%nsym>1) nkcalc=Kmesh%nibz 2141 nkcalc=Kmesh%nbz 2142 !if(1==1)then!DEBUG 2143 !open(67,file="test.rhot1",status="REPLACE") 2144 do ikcalc=1,nkcalc ! for the oscillator strengh, spins are identical without SOC 2145 ! if(Sigp%nkptgw/=Kmesh%nbz) then 2146 ! write(msg,'(6a)')ch10,& 2147 !& ' nkptgw and nbz differs: this is not allowed to compute U in cRPA ' 2148 ! ABI_ERROR(msg) 2149 ! endif 2150 write(std_out,*) " ikcalc",ikcalc,size(Sigp%kptgw2bz),nkcalc 2151 2152 if(mod(ikcalc-1,nprocs)==Wfd%my_rank) then 2153 if(cryst%nsym==1) ik_ibz=Kmesh%tab(ikcalc) ! Index of the irred k-point for GW 2154 if(cryst%nsym>1) ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2155 2156 call prep_calc_ucrpa(ik_ibz,ikcalc,itypatcor,ib1,ib2,Cryst,qp_ebands,Sigp,Gsph_x,Vcp,& 2157 & Kmesh,Qmesh,lcor,M1_q_m,Pawtab,Pawang,Paw_pwff,& 2158 & Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,gwx_ngfft,& 2159 & ngfftf,Dtset%prtvol,Dtset%pawcross,Dtset%plowan_compute,rhot1_q_m,wanbz,rhot1) 2160 end if 2161 end do 2162 if (dtset%plowan_compute<10)then 2163 call xmpi_sum(rhot1_q_m,Wfd%comm,ierr) 2164 call xmpi_sum(M1_q_m,Wfd%comm,ierr) 2165 M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol 2166 rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol 2167 ! do pwx=1,sigp%npwx 2168 ! do ibz=1,Qmesh%nibz 2169 ! do ispinor1=1,dtset%nspinor 2170 ! do ispinor2=1,dtset%nspinor 2171 ! do iatom1=1,cryst%nattyp(itypatcor) 2172 ! do im1=1,2*lcor+1 2173 ! do im2=1,2*lcor+1 2174 ! write(67,*)ibz,im1,im2,rhot1_q_m(iatom1,ispinor1,ispinor2,im1,im2,pwx,ibz) 2175 ! enddo 2176 ! enddo 2177 ! enddo 2178 ! enddo 2179 ! enddo 2180 ! enddo 2181 ! enddo 2182 else 2183 call reduce_operwan_realspace(wanbz,rhot1,sigp%npwx,Qmesh%nibz,Wfd%comm,Kmesh%nbz,Wfd%nsppol) 2184 !call cwtime(cpu,wall,gflops,"start") !reduction of rhot1 2185 ! dim=0 2186 ! do pwx=1,sigp%npwx 2187 ! do ibz=1,Qmesh%nibz 2188 ! do spin=1,wanbz%nsppol 2189 ! do ispinor1=1,wanbz%nspinor 2190 ! do ispinor2=1,wanbz%nspinor 2191 ! do iatom1=1,wanbz%natom_wan 2192 ! do iatom2=1,wanbz%natom_wan 2193 ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1) 2194 ! do pos2=1,size(wanbz%nposition(iatom2)%pos,1) 2195 ! do il1=1,wanbz%nbl_atom_wan(iatom1) 2196 ! do il2=1,wanbz%nbl_atom_wan(iatom2) 2197 ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1 2198 ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1 2199 ! dim=dim+1 2200 ! enddo!im2 2201 ! enddo!im1 2202 ! enddo!il2 2203 ! enddo!il1 2204 ! enddo!pos2 2205 ! enddo!pos1 2206 ! enddo!iatom2 2207 ! enddo!iatom1 2208 ! enddo!ispinor2 2209 ! enddo!ispinor1 2210 ! enddo!spin 2211 ! enddo!ibz 2212 ! enddo!pwx 2213 ! ABI_MALLOC(buffer,(dim)) 2214 ! nnn=0 2215 ! do pwx=1,sigp%npwx 2216 ! do ibz=1,Qmesh%nibz 2217 ! do spin=1,wanbz%nsppol 2218 ! do ispinor1=1,wanbz%nspinor 2219 ! do ispinor2=1,wanbz%nspinor 2220 ! do iatom1=1,wanbz%natom_wan 2221 ! do iatom2=1,wanbz%natom_wan 2222 ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1) 2223 ! do pos2=1,size(wanbz%nposition(iatom2)%pos,1) 2224 ! do il1=1,wanbz%nbl_atom_wan(iatom1) 2225 ! do il2=1,wanbz%nbl_atom_wan(iatom2) 2226 ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1 2227 ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1 2228 ! nnn=nnn+1 2229 ! buffer(nnn)=rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2) 2230 ! enddo!im2 2231 ! enddo!im1 2232 ! enddo!il2 2233 ! enddo!il1 2234 ! enddo!pos2 2235 ! enddo!pos1 2236 ! enddo!iatom2 2237 ! enddo!iatom1 2238 ! enddo!ispinor2 2239 ! enddo!ispinor1 2240 ! enddo!spin 2241 ! enddo!ibz 2242 ! enddo!pwx 2243 ! call xmpi_barrier(Wfd%comm) 2244 ! call xmpi_sum(buffer,Wfd%comm,ierr) 2245 ! call xmpi_barrier(Wfd%comm) 2246 ! buffer=buffer/Kmesh%nbz/Wfd%nsppol 2247 ! nnn=0 2248 ! do pwx=1,sigp%npwx 2249 ! do ibz=1,Qmesh%nibz 2250 ! do spin=1,wanbz%nsppol 2251 ! do ispinor1=1,wanbz%nspinor 2252 ! do ispinor2=1,wanbz%nspinor 2253 ! do iatom1=1,wanbz%natom_wan 2254 ! do iatom2=1,wanbz%natom_wan 2255 ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1) 2256 ! do pos2=1,size(wanbz%nposition(iatom2)%pos,1) 2257 ! do il1=1,wanbz%nbl_atom_wan(iatom1) 2258 ! do il2=1,wanbz%nbl_atom_wan(iatom2) 2259 ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1 2260 ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1 2261 ! nnn=nnn+1 2262 ! rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)=buffer(nnn) 2263 ! !write(67,*)ibz,im1,im2,rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2) 2264 ! enddo!im2 2265 ! enddo!im1 2266 ! enddo!il2 2267 ! enddo!il1 2268 ! enddo!pos2 2269 ! enddo!pos1 2270 ! enddo!iatom2 2271 ! enddo!iatom1 2272 ! enddo!ispinor2 2273 ! enddo!ispinor1 2274 ! enddo!spin 2275 ! enddo!ibz 2276 ! enddo!pwx 2277 ! ABI_FREE(buffer) 2278 endif 2279 2280 !close(67) 2281 ! call xmpi_barrier(Wfd%comm) 2282 ! call cwtime(cpu,wall,gflops,"stop")!reduction of rhot1 2283 ! write(6,*)cpu,wall,gflops 2284 ! if(Cryst%nsym==1) then 2285 ! M1_q_m=M1_q_m/Kmesh%nbz/Wfd%nsppol 2286 ! rhot1_q_m=rhot1_q_m/Kmesh%nbz/Wfd%nsppol 2287 ! endif 2288 ! Calculation of U in cRPA: need to treat with a different cutoff the 2289 ! bare coulomb and the screened coulomb interaction. 2290 ! else !DEBUG 2291 ! write(6,*)"DEBUGGING calc_ucrpa" 2292 ! open(67,file="test.rhot1",status="OLD") 2293 ! rewind(67) 2294 ! do pwx=1,sigp%npwx 2295 ! do ibz=1,Qmesh%nibz 2296 ! do spin=1,wanbz%nsppol 2297 ! do ispinor1=1,wanbz%nspinor 2298 ! do ispinor2=1,wanbz%nspinor 2299 ! do iatom1=1,wanbz%natom_wan 2300 ! do iatom2=1,wanbz%natom_wan 2301 ! do pos1=1,size(wanbz%nposition(iatom1)%pos,1) 2302 ! do pos2=1,size(wanbz%nposition(iatom2)%pos,1) 2303 ! do il1=1,wanbz%nbl_atom_wan(iatom1) 2304 ! do il2=1,wanbz%nbl_atom_wan(iatom2) 2305 ! do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1 2306 ! do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1 2307 ! read(67,*)dummy,dummy,dummy,xx 2308 ! rhot1(pwx,ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl(im1,im2,spin,ispinor1,ispinor2)=xx 2309 ! enddo!im2 2310 ! enddo!im1 2311 ! enddo!il2 2312 ! enddo!il1 2313 ! enddo!pos2 2314 ! enddo!pos1 2315 ! enddo!iatom2 2316 ! enddo!iatom1 2317 ! enddo!ispinor2 2318 ! enddo!ispinor1 2319 ! enddo!spin 2320 ! enddo!ibz 2321 ! enddo!pwx 2322 ! close(67) 2323 ! endif!DEBUG 2324 call calc_ucrpa(itypatcor,cryst,Kmesh,lcor,M1_q_m,Qmesh,Er%npwe,sigp%npwx,& 2325 & Cryst%nsym,Sigp%nomegasr,Sigp%minomega_r,Sigp%maxomega_r,ib1,ib2,& 2326 &'Gsum',Cryst%ucvol,Wfd,Er%fname,dtset%plowan_compute,rhot1,wanbz) 2327 2328 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 2329 !Deallocation of wan and rhot1 2330 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 2331 if (dtset%plowan_compute >=10) then 2332 do pwx=1,sigp%npwx 2333 do ibz=1,Qmesh%nibz 2334 call destroy_operwan_realspace(wanbz,rhot1(pwx,ibz)) 2335 end do 2336 end do 2337 ABI_FREE(rhot1) 2338 call destroy_plowannier(wanbz) 2339 endif 2340 ABI_FREE(rhot1_q_m) 2341 ABI_FREE(M1_q_m) 2342 2343 else 2344 ! MRM: sigmak_todo is an array indicating whether the k-point must be computed for GW correction. This array is set to DO IT (to 1) for any GW calculation 2345 ! and it will only change to NOT DO IT (to 0) in case that a density matrix update calc. is used and we are reading checkpoint files 2346 ABI_MALLOC(sigmak_todo,(Wfd%nkibz)) 2347 sigmak_todo(:)=1 2348 2349 ! ====================================================== 2350 ! ==== Initialize arrays for density matrix updates ==== 2351 ! ====================================================== 2352 if (rdm_update) then 2353 ! Note: all subroutines of 70_gw/m_gwrdm.F90 are implemented assuming nsppol == 1 2354 ABI_CHECK(dtset%nsppol == 1, "1-RDM GW correction only implemented for restricted closed-shell calculations!") 2355 2356 ABI_CALLOC(nateigv, (Wfd%mband, Wfd%mband, Wfd%nkibz, Sigp%nsppol)) 2357 ABI_CALLOC(nat_occs, (Wfd%mband, Wfd%nkibz)) 2358 ABI_CALLOC(xrdm_k_full, (b1gw:b2gw, b1gw:b2gw, Wfd%nkibz)) 2359 2360 write(msg,'(a34,2i9)')' Bands used for the GW 1RDM arrays',b1gw,b2gw 2361 call wrtout(units, msg) 2362 2363 do ik_ibz=1,Wfd%nkibz 2364 do ib=b1gw,b2gw 2365 xrdm_k_full(ib,ib,ik_ibz) = qp_ebands%occ(ib,ik_ibz,1) 2366 end do 2367 do ib=1,Wfd%mband 2368 ! Copy initial occ numbers (in principle 2 or 0 from KS-DFT) 2369 nat_occs(ib,ik_ibz) = qp_ebands%occ(ib,ik_ibz,1) 2370 ! Set to identity matrix 2371 nateigv(ib,ib,ik_ibz,1) = cone 2372 end do 2373 end do 2374 if (readchkprdm) then 2375 gw1rdm_fname = trim(dtfil%fnameabi_chkp_rdm) 2376 call get_chkprdm(Wfd,Kmesh,Sigp,qp_ebands,nat_occs,nateigv,sigmak_todo,my_rank,gw1rdm_fname) 2377 end if 2378 call xmpi_barrier(Wfd%comm) 2379 2380 ! Prepare arrays for the imaginary freq. integration (quadrature) of Sigma_c(iw) 2381 ABI_MALLOC(weights, (Sigp%nomegasi)) 2382 call quadrature_sigma_cw(Sigp,Sr,weights) 2383 end if 2384 2385 ! =================================== 2386 ! ==== Static part (Sigma_x-Vxc) ==== 2387 ! =================================== 2388 do ikcalc=1,Sigp%nkptgw 2389 ! Index of the irred k-point for GW 2390 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 2391 2392 ! Do not compute MELS if the k-point was read from the checkpoint file 2393 ! this IF only affects GW density matrix update! 2394 2395 if (sigmak_todo(ik_ibz) == 1) then 2396 ! min and max band indices for GW corrections (for this k-point) 2397 ib1 = MINVAL(Sigp%minbnd(ikcalc,:)) 2398 ib2 = MAXVAL(Sigp%maxbnd(ikcalc,:)) 2399 call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,qp_ebands,Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,Ltg_k(ikcalc),& 2400 Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,& 2401 gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty) 2402 2403 if (rdm_update) then 2404 ! Only recompute exchange update to the 1-RDM if the point read was broken or not precomputed 2405 2406 ! Compute Sigma_x - Vxc or DELTA Sigma_x - Vxc 2407 ! where DELTA Sigma_x = Sigma_x - hyb_parameter Vx^exact for hyb Functionals. 2408 ! NB: Only restricted closed-shell calcs are implemented here 2409 ABI_CALLOC(pot_k, (ib1:ib2, ib1:ib2)) 2410 ABI_CALLOC(rdm_k, (ib1:ib2, ib1:ib2)) 2411 pot_k(ib1:ib2,ib1:ib2) = Sr%x_mat(ib1:ib2,ib1:ib2,ik_ibz,1) - KS_me%vxcval(ib1:ib2,ib1:ib2,ik_ibz,1) 2412 2413 call calc_rdmx(ib1, ib2, ik_ibz, pot_k, rdm_k, qp_ebands) 2414 2415 ! Update the full 1RDM with the exchange corrected one for this k-point 2416 xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz) = xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz) + rdm_k(ib1:ib2,ib1:ib2) 2417 2418 ! Compute NAT ORBS for exchange corrected 1-RDM 2419 do ib=ib1,ib2 2420 rdm_k(ib,ib) = rdm_k(ib,ib) + qp_ebands%occ(ib,ik_ibz,1) ! Only restricted closed-shell calcs 2421 end do 2422 call natoccs(ib1, ib2, rdm_k, nateigv, nat_occs, qp_ebands, ik_ibz, iinfo=0) ! Only restricted closed-shell calcs 2423 ABI_FREE(pot_k) 2424 ABI_FREE(rdm_k) 2425 end if 2426 else 2427 write(msg,'(a1)') ' ' 2428 call wrtout(std_out,msg) 2429 write(msg,'(a64,i5)') ' Skipping the calc. of Sigma_x - ( Vxc + hyb*Fock) for k-point: ',ik_ibz 2430 call wrtout(std_out,msg) 2431 write(msg,'(a1)') ' ' 2432 call wrtout(std_out,msg) 2433 end if 2434 end do 2435 ! for the time being, do not remove this barrier! 2436 call xmpi_barrier(Wfd%comm) 2437 call timab(421,2,tsec) ! calc_sigx_me 2438 2439 ! ========================================================== 2440 ! ==== Correlation part using the coarse gwc_ngfft mesh ==== 2441 ! ========================================================== 2442 if (mod10/=SIG_HF) then 2443 do ikcalc=1,Sigp%nkptgw 2444 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2445 ib1 = MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2446 ib2 = MAXVAL(Sigp%maxbnd(ikcalc,:)) 2447 2448 ABI_CALLOC(sigcme_k, (nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab)) 2449 2450 if (any(mod10 == [SIG_SEX, SIG_COHSEX])) then 2451 ! Calculate static COHSEX or SEX using the coarse gwc_ngfft mesh. 2452 call cohsex_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Cryst,qp_ebands,Sigp,Sr,Er,Gsph_c,Vcp,Kmesh,Qmesh,& 2453 Ltg_k(ikcalc),Pawtab,Pawang,Paw_pwff,Psps,Wfd,QP_sym,& 2454 gwc_ngfft,Dtset%iomode,Dtset%prtvol,sigcme_k) 2455 else 2456 ! Compute correlated part using the coarse gwc_ngfft mesh. 2457 if (x1rdm/=1 .and. sigmak_todo(ik_ibz)==1) then 2458 ! Do not compute correlation MELS if the k-point was read from the checkpoint file 2459 ! this IF only affects GW density matrix update 2460 call calc_sigc_me(ik_ibz,ikcalc,nomega_sigc,ib1,ib2,Dtset,Cryst,qp_ebands, & 2461 Sigp,Sr,Er,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,& 2462 Ltg_k(ikcalc),PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,& 2463 gwc_ngfft,ngfftf,nfftf,ks_rhor,use_aerhor,ks_aepaw_rhor,sigcme_k) 2464 else 2465 write(msg,'(a1)') ' ' 2466 call wrtout(std_out, msg) 2467 write(msg,'(a44,i5)') ' Skipping the calc. of Sigma_c for k-point: ',ik_ibz 2468 call wrtout(std_out, msg) 2469 write(msg,'(a1)') ' ' 2470 call wrtout(std_out, msg) 2471 end if 2472 end if 2473 2474 ! MRM: compute density matrix numerical correction (freq. integration G0 Sigma_c(iw) G0). 2475 if (rdm_update) then 2476 if (sigmak_todo(ik_ibz) == 1) then 2477 ! Only recompute correlation update to the 1-RDM if the point read was broken or not precomputed 2478 if (x1rdm /= 1) then 2479 ABI_CALLOC(rdm_k, (ib1:ib2,ib1:ib2)) 2480 call calc_rdmc(ib1, ib2, ik_ibz, Sr%omega_i, weights, sigcme_k, qp_ebands, rdm_k) 2481 ! Update the full 1RDM with the GW corrected one for this k-point 2482 ! Only restricted closed-shell calcs 2483 rdm_k(ib1:ib2,ib1:ib2) = xrdm_k_full(ib1:ib2,ib1:ib2,ik_ibz) + rdm_k(ib1:ib2,ib1:ib2) 2484 ! Compute nat orbs and occ numbers at k-point ik_ibz 2485 call natoccs(ib1, ib2, rdm_k, nateigv, nat_occs, qp_ebands, ik_ibz, iinfo=1) 2486 ABI_FREE(rdm_k) 2487 endif 2488 ! Print the checkpoint file if required 2489 if (prtchkprdm) then 2490 gw1rdm_fname = trim(dtfil%fnameabo_chkp_rdm) 2491 call print_chkprdm(Wfd,nat_occs,nateigv,ik_ibz,my_rank,gw1rdm_fname) 2492 end if 2493 end if 2494 else 2495 sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) = sigcme_k 2496 end if 2497 ABI_FREE(sigcme_k) 2498 end do 2499 call xmpi_barrier(Wfd%comm) 2500 if (rdm_update) then 2501 ABI_FREE(xrdm_k_full) 2502 ABI_FREE(weights) 2503 end if 2504 end if 2505 2506 call xmpi_barrier(Wfd%comm) 2507 ABI_FREE(sigmak_todo) 2508 2509 ! 1) Print WFK and DEN files. 2510 ! 2) Build band corrections and compute new energies. 2511 if (rdm_update) then 2512 ABI_CALLOC(gw_rhor, (nfftf, dtset%nspden)) 2513 ! 2514 ! NRM WARNING: only the master has bands on Wfd_nato_master so it prints everything and computes gw_rhor 2515 ! 2516 ! All procs. update the qp_ebands and the Hdr_sigma 2517 call update_hdr_bst(Wfd_nato_master, nat_occs, b1gw, b2gw, qp_ebands, Hdr_sigma, Dtset%ngfft(1:3)) 2518 2519 ! Compute unit cell (averaged) occ = \sum _k weight_k occ_k 2520 call print_tot_occ(qp_ebands) 2521 2522 if (my_rank == master) then 2523 call Wfd_nato_master%rotate(cryst, nateigv, bmask=bdm_mask) ! Let it use bdm_mask and build NOs 2524 call Wfd_nato_master%mkrho(cryst, psps, qp_ebands, ngfftf, nfftf, gw_rhor) ! Construct the density 2525 if (dtset%prtwf == 1) then 2526 ! Print WFK file, here qp_ebands contains nat. orb. occs. 2527 call Wfd_nato_master%write_wfk(Hdr_sigma, qp_ebands, dtfil%fnameabo_wfk, wfknocheck=.True.) 2528 end if 2529 if (dtset%prtden == 1) then 2530 ! Print DEN file 2531 call fftdatar_write("density",dtfil%fnameabo_den,dtset%iomode,Hdr_sigma,& 2532 Cryst,ngfftf,cplex1,nfftf,dtset%nspden,gw_rhor,mpi_enreg_seq,ebands=qp_ebands) 2533 end if 2534 end if 2535 call xmpi_bcast(gw_rhor,master, Wfd%comm, ierr) 2536 2537 ! We no longer need Wfd_nato_master. 2538 call Wfd_nato_master%free() 2539 call xmpi_barrier(Wfd%comm) 2540 2541 ABI_FREE(bdm_mask) ! The master already used bdm_mask 2542 ABI_FREE(nat_occs) ! Occs were already placed in qp_ebands 2543 call em1results_free(Er) ! We no longer need Er for GW@KS-DFT 1RDM but we may need space on the RAM memory 2544 2545 if (gw1rdm==2 .and. Sigp%nkptgw==Wfd%nkibz) then 2546 ! Compute energies only if all k-points are available 2547 ! We need the hole 1-RDM to build Fock[GW.1RDM]! 2548 ABI_CALLOC(old_ks_purex, (b1gw:b2gw, Sigp%nkptgw)) 2549 ABI_CALLOC(new_hartr, (b1gw:b2gw, Sigp%nkptgw)) 2550 ABI_CALLOC(gw_rhog, (2, nfftf)) 2551 ABI_CALLOC(gw_vhartr, (nfftf)) 2552 ! 2553 ! A) Compute Evext = int rho(r) vext(r) dr -> simply dot product on the FFT grid 2554 ! Only restricted closed-shell calcs 2555 ! 2556 den_int = sum(gw_rhor(:,1)) * ucvol / nfftf 2557 evext_energy = sum(gw_rhor(:,1) * vpsp(:)) * ucvol / nfftf 2558 ! 2559 ! B) Coulomb <KS_i|Vh[NO]|KS_j> 2560 ! 2561 ! FFT to build gw_rhog 2562 call fourdp(1, gw_rhog, gw_rhor(:,1), -1, MPI_enreg_seq, nfftf, ndat1, ngfftf, tim_fourdp5) 2563 ecutf = dtset%ecut 2564 if (psps%usepaw == 1) then 2565 ecutf = dtset%pawecutdg 2566 call wrtout(std_out, ch10//' FFT (fine) grid used in PAW GW update:') 2567 end if 2568 call getcut(boxcut, ecutf, gmet, gsqcut, dtset%iboxcut, std_out, k0, ngfftf) 2569 call hartre(1, gsqcut, dtset%icutcoul, psps%usepaw, MPI_enreg_seq, nfftf, ngfftf, dtset%nkpt, dtset%rcut, & 2570 gw_rhog, cryst%rprimd, dtset%vcutgeo, gw_vhartr) 2571 2572 ! Build Vhartree -> gw_vhartr 2573 ABI_ICALLOC(tmp_kstab, (2,Wfd%nkibz,Wfd%nsppol)) 2574 do spin=1,Sigp%nsppol 2575 do ikcalc=1,Sigp%nkptgw 2576 ik_ibz = Kmesh%tab(Sigp%kptgw2bz(ikcalc)) 2577 tmp_kstab(1,ik_ibz,spin)=Sigp%minbnd(ikcalc,spin) 2578 tmp_kstab(2,ik_ibz,spin)=Sigp%maxbnd(ikcalc,spin) 2579 end do 2580 end do 2581 2582 ! Build matrix elements from gw_vhartr -> GW1RDM_me 2583 call calc_vhxc_me(Wfd,KS_mflags,GW1RDM_me,Cryst,Dtset,nfftf,ngfftf,& 2584 ks_vtrial,gw_vhartr,ks_vxc,Psps,Pawtab,KS_paw_an,Pawang,Pawfgrtab,KS_paw_ij,dijexc_core,& 2585 gw_rhor,usexcnhat,ks_nhat,ks_nhatgr,nhatgrdim,tmp_kstab,taur=ks_taur) 2586 2587 call xmpi_barrier(Wfd%comm) 2588 ABI_FREE(tmp_kstab) 2589 ABI_FREE(gw_rhor) 2590 ABI_FREE(gw_rhog) 2591 ABI_FREE(gw_vhartr) 2592 2593 ! Save new <i|Hartree[NO]|i> in KS basis for Delta eik 2594 do ikcalc=1,Sigp%nkptgw 2595 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW 2596 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2597 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2598 do ib=b1gw,b2gw 2599 new_hartr(ib,ikcalc)=GW1RDM_me%vhartree(ib,ib,ik_ibz,1) 2600 end do 2601 end do 2602 ! 2603 ! C) Exchange <KS_i|hyb*K^RANGE-SEP?[KS]|KS_j> and save it in old_ks_purex 2604 ! 2605 ! Build Vcp_ks and Vcp_full 2606 call setup_vcp(Vcp_ks,Vcp_full,Dtset,Gsph_x,Gsph_c,Cryst,Qmesh,Kmesh,coef_hyb,comm) 2607 call xmpi_barrier(Wfd%comm) 2608 if (coef_hyb > tol8) then 2609 do ikcalc=1,Sigp%nkptgw 2610 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW 2611 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2612 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2613 ! Notice we need KS occs and band energy diffs 2614 call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,ks_ebands,Sigp,Sr,Gsph_x,Vcp_ks,Kmesh,Qmesh,Ltg_k(ikcalc),& 2615 Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,QP_sym,& 2616 gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty) 2617 2618 ! Build <KS_i|RS?_Hyb?_Sigma_x[KS]|KS_j> matrix (Use Wfd) 2619 call xmpi_barrier(Wfd%comm) 2620 do ib=b1gw,b2gw 2621 ! Save old alpha*<i|K[KS]|i> from the GS calc. for Delta eik 2622 old_ks_purex(ib,ikcalc)=Sr%x_mat(ib,ib,ik_ibz,1) 2623 end do 2624 end do 2625 endif 2626 call xmpi_barrier(Wfd%comm) 2627 call Vcp_ks%free() 2628 2629 ! Use the Wfd with 'new name' (Wfd_nato_all) because it will contain the GW 1RDM nat. orbs. (bands) 2630 ABI_COMMENT("From now on, the Wfd bands will contain the nat. orbs. ones") 2631 Wfd_nato_all => Wfd 2632 ! Let rotate build the NOs in Wfd_nato_all (KS->NO) 2633 call Wfd_nato_all%rotate(Cryst, nateigv) 2634 call xmpi_barrier(Wfd%comm) 2635 ! 2636 ! D) Non-local <NO_i|Vnl|NO_i> terms [saved on nl_bks(band,k,spin)] 2637 ! 2638 call Wfd_nato_all%get_nl_me(Cryst,Psps,Pawtab,bdm2_mask,nl_bks) 2639 evextnl_energy=zero 2640 do spin=1,Sigp%nsppol 2641 do ikcalc=1,Sigp%nkptgw 2642 do ib=Sr%b1gw,Sr%b2gw 2643 evextnl_energy=evextnl_energy+Kmesh%wt(ikcalc)*qp_ebands%occ(ib,ikcalc,spin)*nl_bks(ib,ikcalc,spin) 2644 end do 2645 end do 2646 end do 2647 call xmpi_barrier(Wfd%comm) 2648 ABI_FREE(nl_bks) 2649 ABI_FREE(bdm2_mask) 2650 ! 2651 ! E) Exchange <NO_i|K[NO]|NO_j> 2652 ! 2653 ! Get the matrix elements and save them on Sr%x_mat 2654 tol_empty=0.0001 ! Allow lower occ numbers 2655 do ikcalc=1,Sigp%nkptgw 2656 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW 2657 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2658 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2659 2660 ! Build <NO_i|Sigma_x[NO]|NO_j> matrix 2661 call calc_sigx_me(ik_ibz,ikcalc,ib1,ib2,Cryst,qp_ebands,Sigp,Sr,Gsph_x,Vcp_full,Kmesh,Qmesh,Ltg_k(ikcalc),& 2662 Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd_nato_all,Wfdf,QP_sym,& 2663 gwx_ngfft,ngfftf,Dtset%prtvol,Dtset%pawcross,tol_empty) 2664 end do 2665 tol_empty=0.01 ! Recover standard value for tolerance on the occ numbers 2666 call xmpi_barrier(Wfd%comm) 2667 call Vcp_full%free() 2668 ! Save the new total exchange energy Ex = Ex[GW.1RDM] 2669 ex_energy=sigma_get_exene(Sr,Kmesh,qp_ebands) 2670 ! Save the new total exchange-correlation MBB energy Exc = Exc^MBB[GW.1RDM] 2671 exc_mbb_energy=sigma_get_excene(Sr,Kmesh,qp_ebands) 2672 2673 ! Transform <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j>, 2674 ! <KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j>, 2675 ! and <KS_i|T|KS_j> -> <NO_i|T|NO_j> 2676 ! 2677 call change_matrix(Sigp,Sr,GW1RDM_me,Kmesh,nateigv) 2678 call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data 2679 2680 ! Print Delta eik and band (state) energies 2681 ! 2682 call print_band_energies(b1gw,b2gw,Sr,Sigp,KS_me,Kmesh,ks_ebands,new_hartr,old_ks_purex) 2683 call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data 2684 ! 2685 ! Print the updated total energy and all energy components 2686 ! 2687 eh_energy=mels_get_haene(Sr,GW1RDM_me,Kmesh,qp_ebands) 2688 ekin_energy=mels_get_kiene(Sr,GW1RDM_me,Kmesh,qp_ebands) 2689 ! SD 2-RDM 2690 etot_sd=ekin_energy+evext_energy+evextnl_energy+QP_energies%e_corepsp+QP_energies%e_ewald+eh_energy+ex_energy 2691 ! MBB 2-RDM 2692 etot_mbb=ekin_energy+evext_energy+evextnl_energy+QP_energies%e_corepsp+QP_energies%e_ewald+eh_energy+exc_mbb_energy 2693 call print_total_energy(ekin_energy,evext_energy,evextnl_energy,QP_energies%e_corepsp,eh_energy,ex_energy,& 2694 exc_mbb_energy,QP_energies%e_ewald,etot_sd,etot_mbb,den_int) 2695 ! 2696 ! Clean GW1RDM_me and allocated arrays 2697 ! 2698 call xmpi_barrier(Wfd%comm) ! Wait for all Sigma_x to be ready before deallocating data 2699 call GW1RDM_me%free() ! Deallocate GW1RD_me 2700 ABI_FREE(old_ks_purex) 2701 ABI_FREE(new_hartr) 2702 else 2703 ABI_FREE(gw_rhor) 2704 ABI_FREE(bdm2_mask) 2705 end if 2706 ABI_FREE(nateigv) 2707 endif 2708 2709 call xmpi_barrier(Wfd%comm) 2710 2711 ! MRM: skip the rest for rdm_update=true (density matrix update) 2712 if (.not.rdm_update) then 2713 ! ===================================================== 2714 ! ==== Solve Dyson equation storing results in Sr% ==== 2715 ! ===================================================== 2716 ! * Use perturbative approach or AC to find QP corrections. 2717 ! * If qp-GW, diagonalize also H0+Sigma in the KS basis set to get the 2718 ! new QP amplitudes and energies (Sr%eigvec_qp and Sr%en_qp_diago. 2719 ! TODO AC with spinor not implemented yet. 2720 ! TODO Diagonalization of Sigma+hhartre with AC is wrong. 2721 ! 2722 call timab(425,1,tsec) ! solve_dyson 2723 do ikcalc=1,Sigp%nkptgw 2724 ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irred k-point for GW 2725 ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point) 2726 ib2=MAXVAL(Sigp%maxbnd(ikcalc,:)) 2727 2728 ABI_MALLOC(sigcme_k,(nomega_sigc,ib2-ib1+1,ib2-ib1+1,Sigp%nsppol*Sigp%nsig_ab)) 2729 sigcme_k=sigcme(:,ib1:ib2,ib1:ib2,ikcalc,:) 2730 2731 call solve_dyson(ikcalc,ib1,ib2,nomega_sigc,Sigp,Kmesh,sigcme_k,qp_ebands%eig,Sr,Dtset%prtvol,Dtfil,Wfd%comm) 2732 ABI_FREE(sigcme_k) 2733 ! 2734 ! Calculate direct gap for each spin and print out final results. 2735 ! We use the valence index of the KS system because we still do not know 2736 ! all the QP corrections. Ideally one should use the QP valence index 2737 do spin=1,Sigp%nsppol 2738 if (Sigp%maxbnd(ikcalc,spin) >= ks_vbik(ik_ibz,spin)+1 .and. & 2739 Sigp%minbnd(ikcalc,spin) <= ks_vbik(ik_ibz,spin) ) then 2740 ks_iv=ks_vbik(ik_ibz,spin) 2741 Sr%egwgap (ik_ibz,spin)= Sr%egw(ks_iv+1,ik_ibz,spin) - Sr%egw(ks_iv,ik_ibz,spin) 2742 Sr%degwgap(ik_ibz,spin)= Sr%degw(ks_iv+1,ik_ibz,spin) - Sr%degw(ks_iv,ik_ibz,spin) 2743 else 2744 ! The "gap" cannot be computed 2745 Sr%e0gap(ik_ibz,spin)=zero; Sr%egwgap (ik_ibz,spin)=zero; Sr%degwgap(ik_ibz,spin)=zero 2746 end if 2747 end do 2748 2749 if (wfd%my_rank == master) call write_sigma_results(ikcalc,ik_ibz,Sigp,Sr,ks_ebands) 2750 end do !ikcalc 2751 2752 call timab(425,2,tsec) ! solve_dyson 2753 call timab(426,1,tsec) ! finalize 2754 2755 ! Update the energies in qp_ebands 2756 ! If QPSCGW, use diagonalized eigenvalues otherwise perturbative results. 2757 if (gwcalctyp>=10) then 2758 do ib=1,Sigp%nbnds 2759 qp_ebands%eig(ib,:,:)=Sr%en_qp_diago(ib,:,:) 2760 end do 2761 else 2762 qp_ebands%eig=Sr%egw 2763 end if 2764 2765 ! ================================================================================ 2766 ! ==== This part is done only if all k-points in the IBZ have been calculated ==== 2767 ! ================================================================================ 2768 if (Sigp%nkptgw==Kmesh%nibz) then 2769 2770 ! Recalculate new occupations and Fermi level. 2771 call ebands_update_occ(qp_ebands,Dtset%spinmagntarget,prtvol=Dtset%prtvol) 2772 qp_vbik(:,:) = ebands_get_valence_idx(qp_ebands) 2773 2774 write(msg,'(2a,3x,2(es16.6,a))')ch10,' New Fermi energy : ',qp_ebands%fermie,' Ha ,',qp_ebands%fermie*Ha_eV,' eV' 2775 call wrtout(units, msg) 2776 2777 ! === If all k-points and all occupied bands are calculated, output EXX === 2778 ! FIXME here be careful about the check on ks_vbik in case of metals 2779 ! if (my_rank==master.and.Sigp%nkptgw==Kmesh%nibz.and.ALL(Sigp%minbnd(:)==1).and.ALL(Sigp%maxbnd(:)>=MAXVAL(nbv(:)))) then 2780 ! if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=MAXVAL(MAXVAL(ks_vbik(:,:),DIM=1))) ) then 2781 if (ALL(Sigp%minbnd==1).and. ALL(Sigp%maxbnd>=ks_vbik) ) then ! FIXME here the two arrays use a different indexing. 2782 2783 ex_energy = sigma_get_exene(Sr,Kmesh,qp_ebands) 2784 write(msg,'(a,2(es16.6,a))')' New Exchange energy : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV' 2785 call wrtout(units, msg) 2786 end if 2787 2788 ! Report the QP gaps (Fundamental and direct) 2789 call ebands_report_gap(qp_ebands,header='QP Band Gaps',unit=ab_out) 2790 2791 ! Band structure interpolation from QP energies computed on the k-mesh. 2792 if (nint(dtset%einterp(1)) /= 0 .and. all(sigp%minbdgw == sigp%minbnd) .and. all(sigp%maxbdgw == sigp%maxbnd)) then 2793 call ebands_interpolate_kpath(qp_ebands, dtset, cryst, [sigp%minbdgw, sigp%maxbdgw], dtfil%filnam_ds(4), comm) 2794 end if 2795 end if ! Sigp%nkptgw==Kmesh%nibz 2796 ! 2797 ! Write SCF data in case of self-consistent calculation === 2798 ! Save Sr%en_qp_diago, Sr%eigvec_qp and m_ks_to_qp in the _QPS file. 2799 ! Note that in the first iteration qp_rhor contains KS rhor, then the mixed rhor. 2800 if (gwcalctyp>=10) then 2801 ! Calculate the new m_ks_to_qp 2802 call updt_m_ks_to_qp(Sigp,Kmesh,nscf,Sr,Sr%m_ks_to_qp) 2803 2804 if (wfd%my_rank == master) then 2805 call wrqps(Dtfil%fnameabo_qps,Sigp,Cryst,Kmesh,Psps,Pawtab,QP_Pawrhoij,& 2806 Dtset%nspden,nscf,nfftf,ngfftf,Sr,qp_ebands,Sr%m_ks_to_qp,qp_rhor) 2807 end if 2808 2809 ! Report the MAX variation for each kptgw and spin 2810 call wrtout(ab_out,ch10//' Convergence of QP corrections ') 2811 do spin=1,Sigp%nsppol 2812 write(msg,'(a,i2,a)')' >>>>> For spin ',spin,' <<<<< ' 2813 call wrtout(ab_out, msg) 2814 do ikcalc=1,Sigp%nkptgw 2815 ib1 = Sigp%minbnd(ikcalc,spin); ib2 = Sigp%maxbnd(ikcalc,spin) 2816 ik_bz = Sigp%kptgw2bz(ikcalc); ik_ibz = Kmesh%tab(ik_bz) 2817 ii = imax_loc(ABS(Sr%degw(ib1:ib2,ik_ibz,spin))) 2818 max_degw = Sr%degw(ii,ik_ibz,spin) 2819 write(msg,('(a,i3,a,2f8.3,a,i3)'))'. kptgw no:',ikcalc,'; Maximum DeltaE = (',max_degw*Ha_eV,') for band index:',ii 2820 call wrtout(ab_out, msg) 2821 end do 2822 end do 2823 end if 2824 2825 ! ============================================ 2826 ! ==== Save the GW results in NETCDF file ==== 2827 ! ============================================ 2828 if (wfd%my_rank == master) then 2829 NCF_CHECK(nctk_open_create(ncid, strcat(dtfil%filnam_ds(4), '_SIGRES.nc'), xmpi_comm_self)) 2830 NCF_CHECK(nctk_defnwrite_ivars(ncid, ["sigres_version"], [1])) 2831 NCF_CHECK(cryst%ncwrite(ncid)) 2832 NCF_CHECK(ebands_ncwrite(ks_ebands, ncid)) 2833 NCF_CHECK(sigma_ncwrite(Sigp, Er, Sr, ncid)) ! WARNING!! If gw1rdm>0 then Er is no longer present!! 2834 ! Add qp_rhor. Note that qp_rhor == ks_rhor if wavefunctions are not updated. 2835 !ncerr = nctk_write_datar("qp_rhor",path,ngfft,cplex,nfft,nspden,& 2836 ! comm_fft,fftn3_distrib,ffti3_local,datar,action) 2837 NCF_CHECK(nf90_close(ncid)) 2838 end if 2839 end if ! MRM skipped if GW density matrix update 2840 end if ! ucrpa 2841 2842 !----------------------------- END OF THE CALCULATION ------------------------ 2843 ! 2844 !===================== 2845 !==== Close Files ==== 2846 !===================== 2847 if (wfd%my_rank == master) then 2848 close(unt_gw ) 2849 close(unt_gwdiag) 2850 close(unt_sig) 2851 close(unt_sgr) 2852 close(unt_sigc) 2853 if (mod10==SIG_GW_AC) close(unt_sgm) 2854 end if 2855 ! 2856 !=============================================== 2857 !==== Free arrays and local data structures ==== 2858 !=============================================== 2859 ABI_FREE(ks_vbik) 2860 ABI_FREE(qp_vbik) 2861 ABI_FREE(ph1d) 2862 ABI_FREE(ph1df) 2863 ABI_FREE(qp_rhor) 2864 ABI_FREE(ks_rhor) 2865 ABI_FREE(ks_rhog) 2866 ABI_FREE(qp_taur) 2867 ABI_FREE(ks_taur) 2868 ABI_FREE(ks_vhartr) 2869 ABI_FREE(ks_vtrial) 2870 ABI_FREE(vpsp) 2871 ABI_FREE(ks_vxc) 2872 ABI_FREE(xccc3d) 2873 ABI_FREE(grchempottn) 2874 ABI_FREE(grewtn) 2875 ABI_FREE(grvdw) 2876 if(.not.rdm_update) then 2877 ABI_FREE(sigcme) 2878 end if 2879 2880 ABI_SFREE(kxc) 2881 ABI_SFREE(qp_vtrial) 2882 ABI_FREE(ks_nhat) 2883 ABI_FREE(ks_nhatgr) 2884 ABI_FREE(dijexc_core) 2885 ABI_FREE(ks_aepaw_rhor) 2886 call pawfgr_destroy(Pawfgr) 2887 2888 ! Deallocation for PAW. 2889 if (Dtset%usepaw==1) then 2890 call pawrhoij_free(KS_Pawrhoij) 2891 ABI_FREE(KS_Pawrhoij) 2892 call pawfgrtab_free(Pawfgrtab) 2893 call paw_ij_free(KS_paw_ij) 2894 call paw_an_free(KS_paw_an) 2895 call pawpwff_free(Paw_pwff) 2896 if (gwcalctyp>=10) then 2897 call pawrhoij_free(QP_pawrhoij) 2898 call paw_ij_free(QP_paw_ij) 2899 call paw_an_free(QP_paw_an) 2900 end if 2901 if (Dtset%pawcross==1) then 2902 call paw_pwaves_lmn_free(Paw_onsite) 2903 call wfdf%free() 2904 end if 2905 end if 2906 ABI_FREE(Paw_onsite) 2907 ABI_FREE(Pawfgrtab) 2908 ABI_FREE(Paw_pwff) 2909 ABI_FREE(KS_paw_ij) 2910 ABI_FREE(KS_paw_an) 2911 if (gwcalctyp>=10) then 2912 ABI_FREE(QP_pawrhoij) 2913 ABI_FREE(QP_paw_an) 2914 ABI_FREE(QP_paw_ij) 2915 end if 2916 2917 call wfd%free() 2918 call destroy_mpi_enreg(MPI_enreg_seq) 2919 call littlegroup_free(Ltg_k) 2920 ABI_FREE(Ltg_k) 2921 call Kmesh%free() 2922 call Qmesh%free() 2923 call Gsph_Max%free() 2924 call Gsph_x%free() 2925 call Gsph_c%free() 2926 call Vcp%free() 2927 call cryst%free() 2928 call sigma_free(Sr) 2929 if(.not.rdm_update) then 2930 call em1results_free(Er) 2931 end if 2932 call ppm_free(PPm) 2933 call Hdr_sigma%free() 2934 call Hdr_wfk%free() 2935 call ebands_free(ks_ebands) 2936 call ebands_free(qp_ebands) 2937 call KS_me%free() 2938 call esymm_free(KS_sym) 2939 ABI_FREE(KS_sym) 2940 2941 if (Sigp%symsigma==1.and.gwcalctyp>=20) then 2942 call esymm_free(QP_sym) 2943 ABI_FREE(QP_sym) 2944 end if 2945 2946 call Sigp%free() 2947 2948 call timab(426,2,tsec) ! finalize 2949 call timab(401,2,tsec) 2950 2951 DBG_EXIT('COLL') 2952 2953 end subroutine sigma