TABLE OF CONTENTS
ABINIT/dfpt_wfkfermi [ Functions ]
NAME
dfpt_wfkfermi
FUNCTION
This routine computes the partial Fermi-level density at a given k-point, and the fixed contribution to the 1st-order Fermi energy (nonlocal and kinetic)
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (DRH, XG, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
cplex=1 if rhoaug is real, 2 if rhoaug is complex
cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
projected with non-local projectors: cprj=<p_i|Cnk>
cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q
projected with non-local projectors: cprjq=<p_i|Cnk+q>
dtfil <type(datafiles_type)>=variables related to files
eig0_k(nband_k)=GS eigenvalues at k (hartree)
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
ibg=shift to be applied on the location of data in the array cprj
ibgq=shift to be applied on the location of data in the array cprjq
icg=shift to be applied on the location of data in the array cg
icgq=shift to be applied on the location of data in the array cgq
idir=direction of the current perturbation
ikpt=number of the k-point
ipert=type of the perturbation
isppol=1 for unpolarized, 2 for spin-polarized
kptopt=option for the generation of k points
mband=maximum number of bands
mcgq=second dimension of the cgq array
mcprjq=second dimension of the cprjq array
mkmem =number of k points treated by this node.
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw or wfs at k
nband_k=number of bands at this k point for that spin polarization
ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
prtvol=control print volume and debugging output
rf_hamkq <type(gs_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
rhoaug(cplex*n4,n5,n6)= density in electrons/bohr**3,
on the augmented fft grid. (cumulative, so input as well as output)
rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
if this ratio has been attributed to the band n (second argument), zero otherwise
wtk_k=weight assigned to the k point.
OUTPUT
eig1_k(2*nband_k**2)=first-order eigenvalues (hartree) fe1fixed_k(nband_k)=contribution to 1st-order Fermi energy from changes of occupation from all bands at this k point. fe1norm_k(nband_k)=contribution to normalization for above rhoaug(cplex*n4,n5,n6)= Fermi-level density in electrons/bohr**3, on the augmented fft grid. (cumulative, so input as well as output). ==== if (gs_hamkq%usepaw==1) ==== pawrhoijfermi(natom) <type(pawrhoij_type)>= paw rhoij occupancies at Fermi level (cumulative, so input as well as output)
PARENTS
dfpt_rhofermi
CHILDREN
dfpt_accrho,dotprod_g,getgh1c,pawcprj_alloc,pawcprj_axpby,pawcprj_copy pawcprj_free,pawcprj_get,status,timab,wrtout
SOURCE
78 #if defined HAVE_CONFIG_H 79 #include "config.h" 80 #endif 81 82 #include "abi_common.h" 83 84 85 subroutine dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,& 86 & dtfil,eig0_k,eig1_k,fe1fixed_k,fe1norm_k,gs_hamkq,& 87 & ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,& 88 & kptopt,mband,mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,& 89 & npw_k,npw1_k,nspinor,nsppol,occ_k,pawrhoijfermi,prtvol,& 90 & rf_hamkq,rhoaug,rocceig,wtk_k) 91 92 use defs_basis 93 use defs_abitypes 94 use m_profiling_abi 95 use m_errors 96 use m_xmpi 97 use m_cgtools 98 99 use m_pawrhoij, only : pawrhoij_type 100 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_copy, pawcprj_axpby, pawcprj_free 101 use m_hamiltonian, only : gs_hamiltonian_type,rf_hamiltonian_type 102 103 !This section has been created automatically by the script Abilint (TD). 104 !Do not modify the following lines by hand. 105 #undef ABI_FUNC 106 #define ABI_FUNC 'dfpt_wfkfermi' 107 use interfaces_14_hidewrite 108 use interfaces_18_timing 109 use interfaces_32_util 110 use interfaces_66_wfs 111 use interfaces_72_response, except_this_one => dfpt_wfkfermi 112 !End of the abilint section 113 114 implicit none 115 116 !Arguments ------------------------------------ 117 !scalars 118 integer,intent(in) :: cplex,ibg,ibgq,icg,icgq,idir,ikpt 119 integer,intent(in) :: ipert,isppol,kptopt,mband,mcgq,mcprjq,mkmem,mpw,ncpgr 120 integer,intent(in) :: npw1_k,nspinor,nsppol,prtvol 121 integer,intent(inout) :: nband_k,npw_k 122 real(dp),intent(in) :: wtk_k 123 type(MPI_type),intent(in) :: mpi_enreg 124 type(datafiles_type),intent(in) :: dtfil 125 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 126 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq 127 !arrays 128 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),cgq(2,mcgq) 129 real(dp),intent(in) :: eig0_k(nband_k),occ_k(nband_k),rocceig(nband_k,nband_k) 130 real(dp),intent(inout) :: rhoaug(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc) 131 real(dp),intent(inout) :: eig1_k(2*nband_k**2) 132 real(dp),intent(out) :: fe1fixed_k(nband_k) 133 real(dp),intent(out) :: fe1norm_k(nband_k) 134 type(pawcprj_type),intent(in) :: cprj(gs_hamkq%natom,nspinor*mband*mkmem*nsppol*gs_hamkq%usecprj) 135 type(pawcprj_type),intent(in) :: cprjq(gs_hamkq%natom,mcprjq) 136 type(pawrhoij_type),intent(inout) :: pawrhoijfermi(gs_hamkq%natom*gs_hamkq%usepaw) 137 138 !Local variables------------------------------- 139 !scalars 140 integer,parameter :: level=18 141 integer :: berryopt,counter,iband,iexit,ii,indx,iorder_cprj 142 integer :: ipw,me,nkpt_max,optlocal,optnl,opt_accrho,opt_corr 143 integer :: opt_gvnl1,sij_opt,tim_fourwf,tim_getgh1c,usevnl 144 real(dp) :: dotr,lambda,wtband 145 character(len=500) :: msg 146 !arrays 147 real(dp) :: dum_grad_berry(1,1),dum_gvnl1(1,1),dum_gs1(1,1),tsec(2) 148 real(dp),allocatable :: cwave0(:,:),cwaveq(:,:),gh1(:,:) 149 type(pawcprj_type),allocatable :: cwaveprj0(:,:),cwaveprjq(:,:),cwaveprj_tmp(:,:) 150 151 ! ********************************************************************* 152 153 DBG_ENTER('COLL') 154 155 !Check arguments validity 156 if (ipert>gs_hamkq%natom.and.ipert/=gs_hamkq%natom+3.and.ipert/=gs_hamkq%natom+4.and.ipert/=gs_hamkq%natom+5) then !SPr rfmagn deb 157 msg=' wrong ipert argument !' 158 MSG_BUG(msg) 159 end if 160 if (cplex/=1) then 161 MSG_BUG('wrong cplex/=1 argument !') 162 end if 163 164 !Debugging statements 165 call status(0,dtfil%filstat,iexit,level,'enter dfpt_wfkfermi') 166 if(prtvol==-level)then 167 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,'dfpt_wfkfermi : enter' 168 call wrtout(std_out,msg,'PERS') 169 end if 170 nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1 171 172 if(prtvol>2 .or. ikpt<=nkpt_max)then 173 write(msg, '(a,a,i5,2x,a,3f9.5,2x,a)' ) ch10,& 174 & ' Non-SCF iterations; k pt #',ikpt,'k=',gs_hamkq%kpt_k(:),' band residuals:' 175 call wrtout(std_out,msg,'PERS') 176 end if 177 178 !Retrieve parallelism data 179 me=mpi_enreg%me_kpt 180 !Initializations and allocations 181 182 ABI_ALLOCATE(gh1,(2,npw1_k*nspinor)) 183 ABI_ALLOCATE(cwave0,(2,npw_k*nspinor)) 184 ABI_ALLOCATE(cwaveq,(2,npw1_k*nspinor)) 185 iorder_cprj=0 ; eig1_k(:)=zero 186 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 187 ABI_DATATYPE_ALLOCATE(cwaveprj0,(gs_hamkq%natom,nspinor)) 188 ABI_DATATYPE_ALLOCATE(cwaveprjq,(gs_hamkq%natom,nspinor)) 189 call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj) 190 call pawcprj_alloc(cwaveprjq,0,gs_hamkq%dimcprj) 191 else 192 ABI_DATATYPE_ALLOCATE(cwaveprj0,(0,0)) 193 ABI_DATATYPE_ALLOCATE(cwaveprjq,(0,0)) 194 end if 195 !Arguments of getgh1c routine (want only (NL+kin) frozen H(1)) 196 berryopt=0;usevnl=0;sij_opt=-gs_hamkq%usepaw;tim_getgh1c=3 197 optlocal=0;optnl=1;opt_gvnl1=0 198 if(ipert==gs_hamkq%natom+5) optnl=0; ! no 1st order NL in H(1), also no kin, but this will be taken into account later 199 !if(ipert==gs_hamkq%natom+5) optlocal=0; ! 1st order LOCAL potential present 200 201 !Arguments of the dfpt_accrho routine 202 tim_fourwf=5 ; opt_accrho=1 ; opt_corr=0 203 !Null potentially unassigned output variables 204 fe1fixed_k(:)=zero; fe1norm_k(:)=zero 205 206 !Read the npw and kg records of wf files 207 call status(0,dtfil%filstat,iexit,level,'before WffRead') 208 209 call timab(139,1,tsec) 210 211 212 !Loop over bands 213 do iband=1,nband_k 214 counter=100*iband+1 215 216 ! Skip bands not treated by current proc 217 if(mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me) cycle 218 219 if(prtvol>=10)then 220 call status(counter,dtfil%filstat,iexit,level,'loop iband ') 221 end if 222 223 ! Select occupied bands 224 if(abs(occ_k(iband))>tol8.and.abs(rocceig(iband,iband))>tol8)then 225 226 wtband=rocceig(iband,iband)/occ_k(iband) 227 ! Get ground-state wavefunctions at k 228 do ipw=1,npw_k*nspinor 229 cwave0(1,ipw)=cg(1,ipw+(iband-1)*npw_k*nspinor+icg) 230 cwave0(2,ipw)=cg(2,ipw+(iband-1)*npw_k*nspinor+icg) 231 end do 232 233 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 234 ! Read PAW ground state projected WF (cprj) 235 call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,gs_hamkq%natom,iband,ibg,ikpt,iorder_cprj,& 236 & isppol,mband,mkmem,gs_hamkq%natom,1,nband_k,nspinor,nsppol,dtfil%unpaw,& 237 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,& 238 & icpgr=idir,ncpgr=ncpgr) 239 end if 240 241 ! Read ground-state wavefunctions at k+q 242 indx=npw1_k*nspinor*(iband-1)+icgq 243 cwaveq(:,1:npw_k*nspinor)=wtband*cgq(:,1+indx:npw_k*nspinor+indx) 244 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 245 ! Read PAW ground state projected WF (cprj) 246 indx=nspinor*(iband-1)+ibgq 247 call pawcprj_copy(cprjq(:,1+indx:nspinor+indx),cwaveprjq) 248 call pawcprj_axpby(zero,wtband,cwaveprj_tmp,cwaveprjq) 249 end if 250 251 if(prtvol>=10)then 252 call status(0,dtfil%filstat,iexit,level,'after wf read ') 253 end if 254 255 256 ! Apply H^(1)-Esp.S^(1) to Psi^(0) (H(^1)=only (NL+kin) frozen part) 257 lambda=eig0_k(iband) 258 call getgh1c(berryopt,cwave0,cwaveprj0,gh1,dum_grad_berry,dum_gs1,gs_hamkq,dum_gvnl1,& 259 & idir,ipert,lambda,mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamkq,sij_opt,& 260 & tim_getgh1c,usevnl) 261 ! Compute Eig1=<Psi^(0)|H^(1)-Eps.S^(1)|Psi(0)> 262 call dotprod_g(dotr,lambda,gs_hamkq%istwf_k,npw_k*nspinor,1,cwave0,gh1,mpi_enreg%me_g0, & 263 & mpi_enreg%comm_spinorfft) 264 indx=2*iband-1+(iband-1)*2*nband_k 265 eig1_k(indx)=dotr 266 ! Compute the fixed contribution to the 1st-order Fermi energy 267 fe1fixed_k(iband)=two*wtband*eig1_k(indx) 268 fe1norm_k(iband) =two*wtband 269 270 ! Accumulate contribution to density and PAW occupation matrix 271 272 call dfpt_accrho(counter,cplex,cwave0,cwaveq,cwaveq,cwaveprj0,cwaveprjq,dotr,& 273 & dtfil%filstat,gs_hamkq,iband,0,0,isppol,kptopt,mpi_enreg,gs_hamkq%natom,nband_k,ncpgr,& 274 & npw_k,npw1_k,nspinor,occ_k,opt_accrho,pawrhoijfermi,prtvol,rhoaug,tim_fourwf,& 275 & opt_corr,wtk_k) 276 277 end if ! End of non-zero occupation and rocceig 278 279 end do ! End loop over bands 280 281 call timab(139,2,tsec) 282 call timab(130,1,tsec) 283 call status(0,dtfil%filstat,iexit,level,'after loops ') 284 285 ABI_DEALLOCATE(cwave0) 286 ABI_DEALLOCATE(cwaveq) 287 ABI_DEALLOCATE(gh1) 288 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 289 call pawcprj_free(cwaveprj0) 290 call pawcprj_free(cwaveprjq) 291 end if 292 ABI_DATATYPE_DEALLOCATE(cwaveprj0) 293 ABI_DATATYPE_DEALLOCATE(cwaveprjq) 294 295 !Structured debugging : if prtvol=-level, stop here. 296 if(prtvol==-level)then 297 write(msg,'(a,a1,a,i2,a)')' fermie3 : exit prtvol=-',level,', debugging mode => stop ' 298 MSG_ERROR(msg) 299 end if 300 301 call status(0,dtfil%filstat,iexit,level,'exit dfpt_wfkfermi') 302 303 call timab(130,2,tsec) 304 305 DBG_EXIT('COLL') 306 307 end subroutine dfpt_wfkfermi