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