TABLE OF CONTENTS
ABINIT/calc_vhxc_me [ Functions ]
NAME
calc_vhxc_me
FUNCTION
Evaluate the matrix elements of $v_H$ and $v_{xc}$ and $v_U$
both in case of NC pseudopotentials and PAW (DFT+U, presently, is only available in PAW)
The matrix elements of $v_{xc}$ are calculated with and without the core contribution.
The later quantity is required in case of GW calculations.
INPUTS
Wfd <type (wfdgw_t)>=Structure gathering information on the wavefunctions. Mflags: Flags specifying the quantities to be computed. Dtset <type(dataset_type)>=all input variables in this dataset ngfftf(18)contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft nfftf=number of points in the fine FFT mesh (for this processor) Pawtab(Dtset%ntypat*Dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh Pawang <type(pawang_type)>=paw angular mesh and related data Paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid Cryst<crystal_t>=unit cell and symmetries vhartr(nfftf)= Hartree potential in real space on the fine FFT mesh vxc(nfftf,nspden)= xc potential in real space on the fine FFT grid rhor(nfftf,nspden)=density in real space (smooth part if PAW). nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc kstab(2,Wfd%nkibz,Wfd%nsppol)=Table temporary used to be compatible with the old implementation.
OUTPUT
Mels
%kinetic=matrix elements of $t$.
%vxc =matrix elements of $v_{xc}[nv+nc]$.
%vxcval =matrix elements of $v_{xc}[nv]$.
%vxcval_hybrid=matrix elements of $v_{xc}[nv]^{hybrid functional}$.
%vhartr =matrix elements of $v_H$.
%vu =matrix elements of $v_U$.
SIDE EFFECTS
Paw_ij= In case of self-Consistency it is changed. It will contain the new H0 Hamiltonian calculated using the QP densities. The valence contribution to XC is removed.
NOTES
All the quantities ($v_H$, $v_{xc}$ and $\psi$ are evaluated on the "fine" FFT mesh.
In case of calculations with NC pseudopotentials the usual mesh is defined by ecut.
For PAW calculations the dense FFT grid defined by pawecutdg is used
Besides, in case of PAW, the matrix elements of V_hartree do not contain the onsite
contributions due to the Coulomb potential generated by ncore and tncore.
These quantities, as well as the onsite kinetic terms, are stored in Paw_ij%dij0.
SOURCE
113 subroutine calc_vhxc_me(Wfd, Mflags, Mels, Cryst, Dtset, nfftf, ngfftf, & 114 vtrial, vhartr, vxc, Psps, Pawtab, Paw_an, Pawang, Pawfgrtab, Paw_ij, dijexc_core, & 115 rhor, usexcnhat, nhat, nhatgr, nhatgrdim, kstab, & 116 taur) ! optional arguments 117 118 !Arguments ------------------------------------ 119 !scalars 120 integer,intent(in) :: nhatgrdim,usexcnhat,nfftf 121 type(Dataset_type),intent(in) :: Dtset 122 type(Pseudopotential_type),intent(in) :: Psps 123 type(wfdgw_t),target,intent(inout) :: Wfd 124 type(Pawang_type),intent(in) :: Pawang 125 type(crystal_t),intent(in) :: Cryst 126 type(melflags_t),intent(in) :: Mflags 127 type(melements_t),intent(out) :: Mels 128 !arrays 129 integer,intent(in) :: ngfftf(18) 130 integer,intent(in) :: kstab(2, Wfd%nkibz, Wfd%nsppol) 131 real(dp),intent(in) :: vhartr(nfftf), vxc(nfftf, Wfd%nspden), vtrial(nfftf, Wfd%nspden) 132 real(dp),intent(in) :: rhor(nfftf, Wfd%nspden) 133 real(dp),intent(in) :: nhat(nfftf, Wfd%nspden * Wfd%usepaw) 134 real(dp),intent(in) :: nhatgr(nfftf, Wfd%nspden, 3 * nhatgrdim) 135 real(dp),intent(in),optional :: taur(nfftf, Wfd%nspden * Dtset%usekden) 136 real(dp),intent(in) :: dijexc_core(:,:,:) ! (cplex_dij*lmn2_size_max,ndij,Cryst%ntypat) 137 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat * Wfd%usepaw) 138 type(Paw_an_type),intent(in) :: Paw_an(Cryst%natom) 139 type(Paw_ij_type),intent(inout) :: Paw_ij(Cryst%natom) 140 type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom) 141 142 !Local variables------------------------------- 143 !scalars 144 integer :: auxc_ixc,iat,ikc,ik_ibz,ib,jb,is,b_start,b_stop,istwf_k 145 integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,lmn2_size_max 146 integer :: isppol,cplex_dij,npw_k,nspinor,nsppol,nspden,nk_calc,rank 147 integer :: iab,isp1,isp2,ixc_sigma,nsploop,nkxc,option,n3xccc_,nk3xc,my_nbbp,my_nmels 148 real(dp) :: nfftfm1,fact,DijH,enxc_val,enxc_hybrid_val,vxcval_avg,vxcval_hybrid_avg,h0dij,vxc1,vxc1_val,re_p,im_p,dijsigcx 149 complex(dpc) :: cdot 150 logical :: ltest,nmxc 151 character(len=500) :: msg 152 type(MPI_type) :: MPI_enreg_seq 153 type(xcdata_type) :: xcdata,xcdata_hybrid 154 type(wave_t),pointer :: wave_jb, wave_ib 155 !arrays 156 integer,parameter :: spinor_idxs(2,4)=RESHAPE([1,1,2,2,1,2,2,1], [2,4]) 157 integer :: got(Wfd%nproc) 158 integer,allocatable :: kcalc2ibz(:),dimlmn(:),bbp_ks_distrb(:,:,:,:) 159 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 160 real(dp) :: tmp_xc(2,Wfd%nspinor**2),tmp_xcval(2,Wfd%nspinor**2) 161 real(dp) :: tmp_H(2,Wfd%nspinor**2),tmp_U(2,Wfd%nspinor**2) 162 real(dp) :: tmp_h0ij(2,Wfd%nspinor**2),tmp_sigcx(2,Wfd%nspinor**2) 163 real(dp) :: dijU(2),strsxc(6),kpt(3),vxc1ab(2),vxc1ab_val(2) 164 real(dp),allocatable :: kxc_(:,:),xccc3d_(:),vxc_val(:,:),vxc_val_hybrid(:,:) 165 real(dp),allocatable :: kinpw(:),veffh0(:,:) 166 complex(dpc) :: tmp(3) 167 complex(gwpc),ABI_CONTIGUOUS pointer :: ur1_up(:),ur1_dwn(:),ur2_up(:),ur2_dwn(:),cg1(:),cg2(:) 168 complex(gwpc),target,allocatable :: ur1(:),ur2(:) 169 complex(dpc),allocatable :: vxcab(:),vxcab_val(:),vxcab_val_hybrid(:),u1cjg_u2dpc(:),kinwf2(:),veffh0_ab(:) 170 logical,allocatable :: bbp_mask(:,:) 171 type(pawcprj_type),allocatable :: Cprj_b1ks(:,:),Cprj_b2ks(:,:) 172 type(libxc_functional_type) :: xc_funcs_hybrid(2) 173 174 ! ************************************************************************* 175 176 DBG_ENTER("COLL") 177 178 ABI_MALLOC(bbp_mask,(Wfd%mband, Wfd%mband)) 179 180 ! Usually FFT meshes for wavefunctions and potentials are not equal. Two approaches are possible: 181 ! Either we Fourier interpolate potentials on the coarse WF mesh or we FFT the wfs on the dense mesh. 182 ! The later approach is used, more CPU demanding but more accurate. 183 if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst, Psps, ngfftf) 184 185 ! Fake MPI_type for sequential part 186 rank = Wfd%my_rank 187 call initmpi_seq(MPI_enreg_seq) 188 call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',ngfftf(2),ngfftf(3),'all') 189 190 nspinor=Wfd%nspinor; nsppol =Wfd%nsppol; nspden =Wfd%nspden 191 if (nspinor == 2) ABI_WARNING("Remember to ADD SO") 192 193 ! TODO not used for the time being but it should be a standard input of the routine. 194 ! bbks_mask(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol)=Logical mask used to select 195 ! the matrix elements to be calculated. 196 ABI_MALLOC(kcalc2ibz,(Wfd%nkibz)) 197 kcalc2ibz=0 198 199 ! Index in the IBZ of the GW k-points. 200 ! Only these points will be considered. 201 nk_calc=0 202 do ik_ibz=1,Wfd%nkibz 203 if ( ALL(kstab(1,ik_ibz,:)/=0) .and. ALL(kstab(2,ik_ibz,:)/=0) ) then 204 nk_calc=nk_calc+1; kcalc2ibz(nk_calc) = ik_ibz 205 end if 206 end do 207 208 call Mels%init(Mflags, nsppol, nspden, Wfd%nspinor, Wfd%nkibz, Wfd%kibz, kstab) 209 210 if (Mflags%has_lexexch==1) then 211 ABI_ERROR("Local EXX not coded!") 212 end if 213 214 ! Evaluate $v_\xc$ using only the valence charge. 215 call wrtout(std_out," calc_vhxc_braket: calculating v_xc[n_val] (excluding non-linear core corrections)") 216 217 do isppol=1,nsppol 218 write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density rhor = ',MINVAL(rhor(:,isppol)) 219 call wrtout(std_out, msg) 220 if (Wfd%usepaw==1) then 221 write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density nhat = ',MINVAL(nhat(:,isppol)) 222 call wrtout(std_out, msg) 223 write(msg,'(a,i2,a,e16.6)')' For spin ',isppol,' Min density trho-nhat = ',MINVAL(rhor(:,isppol)-nhat(:,isppol)) 224 call wrtout(std_out, msg) 225 write(msg,'(a,i2)')' using usexcnhat = ',usexcnhat 226 call wrtout(std_out, msg) 227 end if 228 end do 229 230 option = 0 ! Only exc, vxc, strsxc 231 nkxc = 0 ! No computation of XC kernel 232 n3xccc_= 0 ! No core 233 nk3xc = 0 ! k3xc not needed 234 nmxc = Dtset%usepaw==1 .and. mod(abs(Dtset%usepawu),10) == 4 235 236 ABI_MALLOC(xccc3d_,(n3xccc_)) 237 ABI_MALLOC(kxc_,(nfftf,nkxc)) 238 ABI_MALLOC(vxc_val,(nfftf,nspden)) 239 240 call xcdata_init(xcdata,dtset=Dtset) 241 242 call rhotoxc(enxc_val,kxc_,MPI_enreg_seq,nfftf,ngfftf,& 243 nhat,Wfd%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc_,option,rhor,Cryst%rprimd,& 244 strsxc,usexcnhat,vxc_val,vxcval_avg,xccc3d_,xcdata,taur=taur) 245 246 ! FABIEN's development 247 ! Hybrid functional treatment 248 if (Mflags%has_vxcval_hybrid ==1 ) then 249 250 call wrtout(std_out,' Hybrid functional xc potential is being set') 251 ixc_sigma=Dtset%ixc_sigma 252 call get_auxc_ixc(auxc_ixc,ixc_sigma) 253 call xcdata_init(xcdata_hybrid,dtset=Dtset,auxc_ixc=auxc_ixc,ixc=ixc_sigma) 254 255 if(ixc_sigma<0)then 256 if(libxc_functionals_check()) then 257 call libxc_functionals_init(ixc_sigma,Dtset%nspden,xc_functionals=xc_funcs_hybrid) 258 ! Do not forget, negative values of hyb_mixing(_sr),hyb_range_* means that they have been user-defined. 259 if (dtset%ixc==-402.or.dtset%ixc==-406.or.dtset%ixc==-427.or.dtset%ixc==-428 .or. dtset%ixc==-456 .or. & 260 & min(Dtset%hyb_mixing,Dtset%hyb_mixing_sr,Dtset%hyb_range_dft,Dtset%hyb_range_fock)<-tol8)then 261 call libxc_functionals_set_hybridparams(hyb_range=abs(Dtset%hyb_range_dft),& 262 & hyb_mixing=abs(Dtset%hyb_mixing),hyb_mixing_sr=abs(Dtset%hyb_mixing_sr),xc_functionals=xc_funcs_hybrid) 263 endif 264 else 265 call wrtout(std_out, 'LIBXC is not present: hybrid functionals are not available') 266 end if 267 end if 268 269 write(msg, '(a, f4.2)') ' Fock fraction = ', max(abs(Dtset%hyb_mixing),abs(Dtset%hyb_mixing_sr)) 270 call wrtout(std_out, msg) 271 write(msg, '(a, f5.2, a)') ' Fock inverse screening length = ',abs(Dtset%hyb_range_dft), ' (bohr^-1)' 272 call wrtout(std_out, msg) 273 274 ABI_MALLOC(vxc_val_hybrid,(nfftf,nspden)) 275 276 if(ixc_sigma<0)then 277 call rhotoxc(enxc_hybrid_val,kxc_,MPI_enreg_seq,nfftf,ngfftf,& 278 nhat,Wfd%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc_,option,rhor,Cryst%rprimd,& 279 strsxc,usexcnhat,vxc_val_hybrid,vxcval_hybrid_avg,xccc3d_,xcdata_hybrid,xc_funcs=xc_funcs_hybrid) 280 call libxc_functionals_end(xc_functionals=xc_funcs_hybrid) 281 else 282 call rhotoxc(enxc_hybrid_val,kxc_,MPI_enreg_seq,nfftf,ngfftf,& 283 nhat,Wfd%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc_,option,rhor,Cryst%rprimd,& 284 strsxc,usexcnhat,vxc_val_hybrid,vxcval_hybrid_avg,xccc3d_,xcdata_hybrid) 285 end if 286 287 endif 288 289 ABI_FREE(xccc3d_) 290 ABI_FREE(kxc_) 291 292 write(msg,'(a,f8.4,2a,f8.4,a)')' E_xc[n_val] = ',enxc_val, ' [Ha]. ','<V_xc[n_val]> = ',vxcval_avg,' [Ha]. ' 293 call wrtout(std_out, msg) 294 295 ! has_hbare uses veffh0. Why only use it with usepaw=1? Let it be available always 296 if (Mflags%has_hbare == 1) then 297 if (Mflags%has_kinetic/=1) then 298 ABI_ERROR("Kinetic energy mels are required for the construction of Hbare mels!") 299 end if 300 ! Effective potential of the bare Hamiltonian: valence term is subtracted. 301 ABI_MALLOC(veffh0,(nfftf,nspden)) 302 veffh0=vtrial-vxc_val 303 !veffh0=vtrial !this is to retrieve the KS Hamiltonian 304 endif 305 306 ! If PAW and qp-SCGW then update Paw_ij and calculate the matrix elements === 307 ! We cannot simply rely on gwcalctyp because I need KS vxc in sigma. 308 if (Wfd%usepaw==1.and.Mflags%has_hbare==1) then 309 ABI_CHECK(Mflags%only_diago==0,"Wrong only_diago") 310 311 call paw_mknewh0(Cryst%natom,nsppol,nspden,nfftf,Dtset%pawspnorb,Dtset%pawprtvol,Cryst,& 312 Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial) 313 314 ! Effective potential of the bare Hamiltonian: valence term is subtracted. 315 veffh0=vtrial-vxc_val 316 !veffh0=vtrial !this is to retrieve the KS Hamiltonian 317 end if 318 319 ! Setup of the hermitian operator vxcab === 320 ! if nspden==4 vxc contains (v^11, v^22, Re[V^12], Im[V^12]. 321 ! Cannot use directly Re and Im since we also need off-diagonal elements. 322 if (wfd%nspden == 4) then 323 ABI_MALLOC(vxcab, (nfftf)) 324 ABI_MALLOC(vxcab_val, (nfftf)) 325 vxcab (:) = DCMPLX(vxc (:,3), vxc (:,4)) 326 vxcab_val(:) = DCMPLX(vxc_val(:,3), vxc_val(:,4)) 327 if (Mflags%has_vxcval_hybrid==1) then 328 ABI_MALLOC(vxcab_val_hybrid,(nfftf)) 329 vxcab_val_hybrid(:)=DCMPLX(vxc_val_hybrid(:,3),vxc_val_hybrid(:,4)) 330 end if 331 if (Mflags%has_hbare==1) then 332 ABI_MALLOC(veffh0_ab,(nfftf)) 333 veffh0_ab(:)=DCMPLX(veffh0(:,3),veffh0(:,4)) 334 end if 335 end if 336 337 ABI_MALLOC(ur1, (nfftf * nspinor)) 338 ABI_MALLOC(ur2, (nfftf * nspinor)) 339 ABI_MALLOC(u1cjg_u2dpc, (nfftf * nspinor)) 340 341 ! Create distribution table for tasks. 342 ! This section is parallelized inside wfd%comm 343 ! as all processors are calling the routine with all GW wavefunctions 344 ! TODO the table can be calculated at each (k,s) to save some memory. 345 got=0; my_nmels=0 346 ABI_MALLOC(bbp_ks_distrb,(Wfd%mband,Wfd%mband,nk_calc,nsppol)) 347 do is=1,nsppol 348 do ikc=1,nk_calc 349 ik_ibz=kcalc2ibz(ikc) 350 bbp_mask=.FALSE. 351 b_start=kstab(1,ik_ibz,is) 352 b_stop =kstab(2,ik_ibz,is) 353 if (Mflags%only_diago==1) then 354 !do jb=b1,b2 355 do jb=b_start,b_stop 356 bbp_mask(jb,jb)=.TRUE. 357 end do 358 else 359 bbp_mask(b_start:b_stop,b_start:b_stop)=.TRUE. 360 end if 361 362 call wfd%distribute_bbp(ik_ibz,is,"Upper",my_nbbp,bbp_ks_distrb(:,:,ikc,is),got,bbp_mask) 363 my_nmels = my_nmels + my_nbbp 364 end do 365 end do 366 ABI_FREE(bbp_mask) 367 368 write(msg,'(a,i0,a)')" Will calculate ",my_nmels," <b,k,s|O|b',k,s> matrix elements in calc_vhxc_me." 369 call wrtout(std_out, msg) 370 371 ! ===================================== 372 ! ==== Loop over required k-points ==== 373 ! ===================================== 374 nfftfm1=one/nfftf 375 376 do is=1,nsppol 377 if (ALL(bbp_ks_distrb(:,:,:,is)/=rank)) CYCLE 378 do ikc=1,nk_calc 379 if (ALL(bbp_ks_distrb(:,:,ikc,is)/=rank)) CYCLE 380 381 ik_ibz = kcalc2ibz(ikc) 382 b_start = kstab(1,ik_ibz,is) 383 b_stop = kstab(2,ik_ibz,is) 384 npw_k = Wfd%Kdata(ik_ibz)%npw 385 kpt = Wfd%kibz(:,ik_ibz) 386 kg_k => Wfd%kdata(ik_ibz)%kg_k 387 istwf_k = wfd%istwfk(ik_ibz) 388 389 ! Calculate |k+G|^2 needed by hbareme and kineticme 390 ! MRM: Solved ecut problem 391 if (Mflags%has_kinetic == 1) then 392 ABI_MALLOC(kinpw, (npw_k)) 393 ABI_MALLOC(kinwf2, (npw_k*nspinor)) 394 call mkkin(Dtset%ecutwfn,Dtset%ecutsm,Dtset%effmass_free,Cryst%gmet,kg_k,kinpw,kpt,npw_k,0,0) 395 where (kinpw>HUGE(zero)*1.d-11) 396 kinpw=zero 397 end where 398 end if 399 400 !do jb=b1,b2 401 do jb=b_start,b_stop 402 if (ALL(bbp_ks_distrb(:,jb,ikc,is)/=rank)) CYCLE 403 404 ABI_CHECK(wfd%get_wave_ptr(jb, ik_ibz, is, wave_jb, msg) == 0, msg) 405 406 if (Mflags%has_kinetic == 1) then 407 cg2 => wave_jb%ug 408 kinwf2(1:npw_k) = cg2(1:npw_k)*kinpw(:) 409 if (nspinor==2) kinwf2(npw_k+1:)=cg2(npw_k+1:)*kinpw(:) 410 end if 411 412 call wfd%get_ur(jb,ik_ibz,is,ur2) 413 414 !do ib=b1,jb ! Upper triangle 415 do ib=b_start,jb 416 if (bbp_ks_distrb(ib,jb,ikc,is)/=rank) CYCLE 417 ! Off-diagonal elements only for QPSCGW. 418 if (Mflags%only_diago==1.and.ib/=jb) CYCLE 419 420 call wfd%get_ur(ib,ik_ibz,is,ur1) 421 u1cjg_u2dpc(:) = CONJG(ur1) *ur2 422 423 if (Mflags%has_vxc == 1) then 424 Mels%vxc(ib, jb, ik_ibz, is) = sum(u1cjg_u2dpc(1:nfftf) * vxc(1:nfftf, is)) * nfftfm1 425 if (wfd%nspinor == 2 .and. wfd%nspden == 1) & 426 Mels%vxc(ib, jb, ik_ibz, 2) = sum(u1cjg_u2dpc(nfftf+1:) * vxc(1:nfftf, is)) * nfftfm1 427 end if 428 if (Mflags%has_vxcval == 1) then 429 Mels%vxcval(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vxc_val(1:nfftf, is)) * nfftfm1 430 if (wfd%nspinor == 2 .and. wfd%nspden == 1) & 431 Mels%vxcval(ib, jb, ik_ibz, 2) = SUM(u1cjg_u2dpc(nfftf+1:) * vxc_val(1:nfftf, is)) * nfftfm1 432 end if 433 if (Mflags%has_vxcval_hybrid == 1) then 434 Mels%vxcval_hybrid(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vxc_val_hybrid(1:nfftf, is)) * nfftfm1 435 if (wfd%nspinor == 2 .and. wfd%nspden == 1) & 436 Mels%vxcval_hybrid(ib, jb, ik_ibz, 2) = SUM(u1cjg_u2dpc(nfftf+1) * vxc_val_hybrid(1:nfftf, is)) * nfftfm1 437 end if 438 if (Mflags%has_vhartree==1) then 439 Mels%vhartree(ib, jb, ik_ibz, is) = SUM(u1cjg_u2dpc(1:nfftf) * vhartr(1:nfftf)) * nfftfm1 440 if (wfd%nspinor == 2 .and. wfd%nspden == 1) & 441 Mels%vhartree(ib, jb, ik_ibz, 2) = SUM(u1cjg_u2dpc(nfftf+1:) * vhartr(1:nfftf)) * nfftfm1 442 end if 443 if (Mflags%has_kinetic==1) then 444 ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, is, wave_ib, msg) == 0, msg) 445 cg1 => wave_ib%ug(1:npw_k) 446 cdot = DOT_PRODUCT(cg1, kinwf2(1:npw_k)) 447 !if (istwf_k /= 1) then 448 ! cdot = two * cdot; if (istwf_k == 2) cdot = cdot - GWPC_CONJG(cg1(1)) * kinwf2(1) 449 !end if 450 Mels%kinetic(ib, jb, ik_ibz, is) = cdot 451 if (wfd%nspinor == 2 .and. wfd%nspden == 1) then 452 cg1 => wave_ib%ug(npw_k+1:) 453 Mels%kinetic(ib, jb, ik_ibz, 2) = DOT_PRODUCT(cg1, kinwf2(npw_k+1:)) 454 end if 455 end if 456 if (Mflags%has_hbare==1) then 457 Mels%hbare(ib, jb, ik_ibz, is) = Mels%kinetic(ib, jb, ik_ibz, is) & 458 + SUM(u1cjg_u2dpc(1:nfftf) * veffh0(1:nfftf, is)) * nfftfm1 459 if (wfd%nspinor == 2 .and. wfd%nspden == 1) then 460 cg1 => wave_ib%ug(npw_k+1:) 461 Mels%hbare(ib, jb, ik_ibz, 2) = & 462 Mels%kinetic(ib, jb, ik_ibz, 2) + SUM(u1cjg_u2dpc(nfftf+1:) * veffh0(1:nfftf, is)) * nfftfm1 463 end if 464 end if 465 466 if (nspinor == 2 .and. wfd%nspden == 4) then 467 ! Here I can skip 21 if ib==jb 468 ur1_up => ur1(1:nfftf) 469 ur1_dwn => ur1(nfftf+1:2*nfftf) 470 ur2_up => ur2(1:nfftf) 471 ur2_dwn => ur2(nfftf+1:2*nfftf) 472 473 if (Mflags%has_kinetic==1) then 474 cg1 => wave_ib%ug(npw_k+1:) 475 tmp(1)=DOT_PRODUCT(cg1,kinwf2(npw_k+1:)) 476 Mels%kinetic(ib,jb,ik_ibz,2 )=tmp(1) 477 Mels%kinetic(ib,jb,ik_ibz,3:4)=czero 478 end if 479 if (Mflags%has_hbare==1) then 480 cg1 => wave_ib%ug(npw_k+1:) 481 tmp(1)=SUM(CONJG(ur1_dwn)*veffh0(:,2)*ur2_dwn)*nfftfm1 + Mels%kinetic(ib,jb,ik_ibz,2) 482 tmp(2)=SUM(CONJG(ur1_dwn)* veffh0_ab(:) *ur2_dwn)*nfftfm1 483 tmp(3)=SUM(CONJG(ur1_dwn)*CONJG(veffh0_ab(:))*ur2_dwn)*nfftfm1 484 Mels%hbare(ib,jb,ik_ibz,2:4)=tmp(:) 485 end if 486 if (Mflags%has_vxc==1) then 487 tmp(1) = SUM(CONJG(ur1_dwn)* vxc(:,2) *ur2_dwn)*nfftfm1 488 tmp(2) = SUM(CONJG(ur1_up )* vxcab(:) *ur2_dwn)*nfftfm1 489 tmp(3) = SUM(CONJG(ur1_dwn)*CONJG(vxcab(:))*ur2_up )*nfftfm1 490 Mels%vxc(ib,jb,ik_ibz,2:4)=tmp(:) 491 end if 492 if (Mflags%has_vxcval==1) then 493 tmp(1) = SUM(CONJG(ur1_dwn)* vxc_val(:,2) *ur2_dwn)*nfftfm1 494 tmp(2) = SUM(CONJG(ur1_up )* vxcab_val(:) *ur2_dwn)*nfftfm1 495 tmp(3) = SUM(CONJG(ur1_dwn)*CONJG(vxcab_val(:))*ur2_up )*nfftfm1 496 Mels%vxcval(ib,jb,ik_ibz,2:4)=tmp(:) 497 end if 498 if (Mflags%has_vxcval_hybrid==1) then 499 tmp(1) = SUM(CONJG(ur1_dwn)* vxc_val_hybrid(:,2) *ur2_dwn)*nfftfm1 500 tmp(2) = SUM(CONJG(ur1_up )* vxcab_val_hybrid(:) *ur2_dwn)*nfftfm1 501 tmp(3) = SUM(CONJG(ur1_dwn)*CONJG(vxcab_val_hybrid(:))*ur2_up )*nfftfm1 502 Mels%vxcval_hybrid(ib,jb,ik_ibz,2:4)=tmp(:) 503 end if 504 if (Mflags%has_vhartree==1) then 505 tmp(1) = SUM(CONJG(ur1_dwn)*vhartr(:)*ur2_dwn)*nfftfm1 506 Mels%vhartree(ib,jb,ik_ibz,2 )=tmp(1) 507 Mels%vhartree(ib,jb,ik_ibz,3:4)=czero 508 end if 509 end if !nspinor==2 510 511 end do !ib 512 end do !jb 513 514 if (Mflags%has_kinetic==1) then 515 ABI_FREE(kinpw) 516 ABI_FREE(kinwf2) 517 end if 518 519 end do !ikc 520 end do !is 521 522 ABI_FREE(ur1) 523 ABI_FREE(ur2) 524 ABI_FREE(vxc_val) 525 ABI_FREE(u1cjg_u2dpc) 526 if(Mflags%has_vxcval_hybrid==1) then 527 ABI_FREE(vxc_val_hybrid) 528 end if 529 if (wfd%nspden == 4) then 530 ABI_FREE(vxcab) 531 ABI_FREE(vxcab_val) 532 if(Mflags%has_vxcval_hybrid==1) then 533 ABI_FREE(vxcab_val_hybrid) 534 end if 535 end if 536 537 if (Mflags%has_hbare==1) then 538 ABI_FREE(veffh0) 539 if (nspinor==2) then 540 ABI_FREE(veffh0_ab) 541 end if 542 end if 543 544 ! ==================================== 545 ! ===== Additional terms for PAW ===== 546 ! ==================================== 547 if (Wfd%usepaw==1) then 548 ! Tests if needed pointers in Paw_ij are allocated. 549 ltest=(allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_hat).and.allocated(Paw_ij(1)%dijxc_val)) 550 ABI_CHECK(ltest,"dijxc, dijxc_hat or dijxc_val not allocated") 551 ABI_CHECK(nspinor == 1, "PAW with nspinor not tested") 552 553 ! For DFT+U 554 do iat=1,Cryst%natom 555 itypat=Cryst%typat(iat) 556 if (Pawtab(itypat)%usepawu/=0) then 557 ltest=(allocated(Paw_ij(iat)%dijU)) 558 ABI_CHECK(ltest,"DFT+U but dijU not allocated") 559 end if 560 end do 561 562 if (Dtset%pawspnorb>0) then 563 ltest=(allocated(Paw_ij(1)%dijso)) 564 ABI_CHECK(ltest,"dijso not allocated") 565 end if 566 567 lmn2_size_max=MAXVAL(Pawtab(:)%lmn2_size) 568 569 if (Mflags%has_sxcore==1) then 570 if ( SIZE(dijexc_core,DIM=1) /= lmn2_size_max & 571 .or. SIZE(dijexc_core,DIM=2) /= 1 & 572 .or. SIZE(dijexc_core,DIM=3) /= Cryst%ntypat ) then 573 ABI_BUG("Wrong sizes in dijexc_core") 574 end if 575 end if 576 577 nsploop=nspinor**2 578 579 ! ==================================== 580 ! === Assemble PAW matrix elements === 581 ! ==================================== 582 ABI_MALLOC(dimlmn, (Cryst%natom)) 583 do iat=1,Cryst%natom 584 dimlmn(iat)=Pawtab(Cryst%typat(iat))%lmn_size 585 end do 586 587 ABI_MALLOC(Cprj_b1ks, (Cryst%natom,nspinor)) 588 ABI_MALLOC(Cprj_b2ks, (Cryst%natom,nspinor)) 589 call pawcprj_alloc(Cprj_b1ks, 0, dimlmn) 590 call pawcprj_alloc(Cprj_b2ks, 0, dimlmn) 591 592 do is=1,nsppol 593 if (ALL(bbp_ks_distrb(:,:,:,is)/=rank)) CYCLE 594 595 ! Loop over required k-points 596 do ikc=1,nk_calc 597 if (ALL(bbp_ks_distrb(:,:,ikc,is)/=rank)) CYCLE 598 ik_ibz=kcalc2ibz(ikc) 599 b_start=kstab(1,ik_ibz,is) 600 b_stop =kstab(2,ik_ibz,is) 601 602 !do jb=b1,b2 603 do jb=b_start,b_stop 604 if (ALL(bbp_ks_distrb(:,jb,ikc,is)/=rank)) CYCLE 605 606 ! Load projected wavefunctions for this k-point, spin and band === 607 ! Cprj are unsorted, full correspondence with xred. See ctocprj.F90!! 608 call wfd%get_cprj(jb,ik_ibz,is,Cryst,Cprj_b2ks,sorted=.FALSE.) 609 610 !do ib=b1,jb ! Upper triangle 611 do ib=b_start,jb 612 if (bbp_ks_distrb(ib,jb,ikc,is)/=rank) CYCLE 613 614 ! Off-diagonal elements only for QPSCGW. 615 if (Mflags%only_diago==1.and.ib/=jb) CYCLE 616 617 call wfd%get_cprj(ib,ik_ibz,is,Cryst,Cprj_b1ks,sorted=.FALSE.) 618 ! 619 ! === Get onsite matrix elements summing over atoms and channels === 620 ! * Spin is external and fixed (1,2) if collinear. 621 ! * if noncollinear loop internally over the four components ab. 622 tmp_xc = zero; tmp_xcval = zero; tmp_H = zero; tmp_U = zero; tmp_h0ij = zero; tmp_sigcx = zero 623 624 do iat=1,Cryst%natom 625 itypat =Cryst%typat(iat) 626 lmn_size =Pawtab(itypat)%lmn_size 627 cplex_dij=Paw_ij(iat)%cplex_dij 628 klmn1=1 629 630 do jlmn=1,lmn_size 631 j0lmn=jlmn*(jlmn-1)/2 632 do ilmn=1,jlmn 633 klmn=j0lmn+ilmn 634 ! TODO Be careful, here I assume that the onsite terms ij are symmetric 635 ! should check the spin-orbit case! 636 fact=one; if (ilmn==jlmn) fact=half 637 638 ! Loop over four components if nspinor==2 639 ! If collinear nsploop==1 640 do iab=1,nsploop 641 isp1=spinor_idxs(1,iab); isp2=spinor_idxs(2,iab) 642 643 re_p= Cprj_b1ks(iat,isp1)%cp(1,ilmn) * Cprj_b2ks(iat,isp2)%cp(1,jlmn) & 644 +Cprj_b1ks(iat,isp1)%cp(2,ilmn) * Cprj_b2ks(iat,isp2)%cp(2,jlmn) & 645 +Cprj_b1ks(iat,isp1)%cp(1,jlmn) * Cprj_b2ks(iat,isp2)%cp(1,ilmn) & 646 +Cprj_b1ks(iat,isp1)%cp(2,jlmn) * Cprj_b2ks(iat,isp2)%cp(2,ilmn) 647 648 im_p= Cprj_b1ks(iat,isp1)%cp(1,ilmn) * Cprj_b2ks(iat,isp2)%cp(2,jlmn) & 649 -Cprj_b1ks(iat,isp1)%cp(2,ilmn) * Cprj_b2ks(iat,isp2)%cp(1,jlmn) & 650 +Cprj_b1ks(iat,isp1)%cp(1,jlmn) * Cprj_b2ks(iat,isp2)%cp(2,ilmn) & 651 -Cprj_b1ks(iat,isp1)%cp(2,jlmn) * Cprj_b2ks(iat,isp2)%cp(1,ilmn) 652 653 ! ================================================== 654 ! === Load onsite matrix elements and accumulate === 655 ! ================================================== 656 if (nspinor==1) then 657 658 if (Mflags%has_hbare==1) then ! * Get new dij of h0 and accumulate. 659 h0dij=Paw_ij(iat)%dij(klmn,is) 660 tmp_h0ij(1,iab)=tmp_h0ij(1,iab) + h0dij*re_p*fact 661 tmp_h0ij(2,iab)=tmp_h0ij(2,iab) + h0dij*im_p*fact 662 end if 663 664 if (Mflags%has_sxcore==1) then ! * Fock operator generated by core electrons. 665 dijsigcx = dijexc_core(klmn,1,itypat) 666 tmp_sigcx(1,iab)=tmp_sigcx(1,iab) + dijsigcx*re_p*fact 667 tmp_sigcx(2,iab)=tmp_sigcx(2,iab) + dijsigcx*im_p*fact 668 end if 669 670 if (Mflags%has_vxc==1) then ! * Accumulate vxc[n1+nc] + vxc[n1+tn+nc]. 671 vxc1 = Paw_ij(iat)%dijxc(klmn,is)+Paw_ij(iat)%dijxc_hat(klmn,is) 672 tmp_xc(1,iab)=tmp_xc(1,iab) + vxc1*re_p*fact 673 tmp_xc(2,iab)=tmp_xc(2,iab) + vxc1*im_p*fact 674 end if 675 676 if (Mflags%has_vxcval==1) then ! * Accumulate valence-only XC. 677 vxc1_val=Paw_ij(iat)%dijxc_val(klmn,is) 678 tmp_xcval(1,1)=tmp_xcval(1,1) + vxc1_val*re_p*fact 679 tmp_xcval(2,1)=tmp_xcval(2,1) + vxc1_val*im_p*fact 680 end if 681 682 if (Mflags%has_vhartree==1) then ! * Accumulate Hartree term of the PAW Hamiltonian. 683 DijH=Paw_ij(iat)%dijhartree(klmn) 684 tmp_H(1,1)=tmp_H(1,1) + DijH*re_p*fact 685 tmp_H(2,1)=tmp_H(2,1) + DijH*im_p*fact 686 end if 687 688 ! Accumulate U term of the PAW Hamiltonian (only onsite AE contribution) 689 if (Mflags%has_vu==1) then 690 if (Pawtab(itypat)%usepawu/=0) then 691 dijU(1)=Paw_ij(iat)%dijU(klmn,is) 692 tmp_U(1,1)=tmp_U(1,1) + dijU(1)*re_p*fact 693 tmp_U(2,1)=tmp_U(2,1) + dijU(1)*im_p*fact 694 end if 695 end if 696 697 else 698 ! Spinorial case 699 700 ! FIXME H0 + spinor not implemented 701 if (Mflags%has_hbare==1.or.Mflags%has_sxcore==1) then 702 ABI_ERROR("not implemented") 703 end if 704 705 if (Mflags%has_vxc==1) then ! * Accumulate vxc[n1+nc] + vxc[n1+tn+nc]. 706 vxc1ab(1) = Paw_ij(iat)%dijxc(klmn1, iab)+Paw_ij(iat)%dijxc_hat(klmn1, iab) 707 vxc1ab(2) = Paw_ij(iat)%dijxc(klmn1+1,iab)+Paw_ij(iat)%dijxc_hat(klmn1+1,iab) 708 tmp_xc(1,iab) = tmp_xc(1,iab) + (vxc1ab(1)*re_p - vxc1ab(2)*im_p)*fact 709 tmp_xc(2,iab) = tmp_xc(2,iab) + (vxc1ab(2)*re_p + vxc1ab(1)*im_p)*fact 710 end if 711 712 if (Mflags%has_vxcval==1) then ! * Accumulate valence-only XC. 713 vxc1ab_val(1) = Paw_ij(iat)%dijxc_val(klmn1, iab) 714 vxc1ab_val(2) = Paw_ij(iat)%dijxc_val(klmn1+1,iab) 715 tmp_xcval(1,iab) = tmp_xcval(1,iab) + (vxc1ab_val(1)*re_p - vxc1ab_val(2)*im_p)*fact 716 tmp_xcval(2,iab) = tmp_xcval(2,iab) + (vxc1ab_val(2)*re_p + vxc1ab_val(1)*im_p)*fact 717 end if 718 719 ! * In GW, dijhartree is always real. 720 if (Mflags%has_vhartree==1) then ! * Accumulate Hartree term of the PAW Hamiltonian. 721 if (iab==1.or.iab==2) then 722 DijH = Paw_ij(iat)%dijhartree(klmn) 723 tmp_H(1,iab) = tmp_H(1,iab) + DijH*re_p*fact 724 tmp_H(2,iab) = tmp_H(2,iab) + DijH*im_p*fact 725 end if 726 end if 727 728 ! TODO "ADD DFT+U and SO" 729 ! check this part 730 if (Mflags%has_vu==1) then 731 if (Pawtab(itypat)%usepawu/=0) then 732 ! Accumulate the U term of the PAW Hamiltonian (only onsite AE contribution) 733 dijU(1)=Paw_ij(iat)%dijU(klmn1 ,iab) 734 dijU(2)=Paw_ij(iat)%dijU(klmn1+1,iab) 735 tmp_U(1,iab) = tmp_U(1,iab) + (dijU(1)*re_p - dijU(2)*im_p)*fact 736 tmp_U(2,iab) = tmp_U(2,iab) + (dijU(2)*re_p + dijU(1)*im_p)*fact 737 end if 738 end if 739 740 end if 741 end do !iab 742 743 klmn1=klmn1+cplex_dij 744 745 end do !ilmn 746 end do !jlmn 747 end do !iat 748 749 ! ======================================== 750 ! ==== Add to plane wave contribution ==== 751 ! ======================================== 752 if (nspinor==1) then 753 754 if (Mflags%has_hbare==1) & 755 & Mels%hbare(ib,jb,ik_ibz,is) = Mels%hbare(ib,jb,ik_ibz,is) + DCMPLX(tmp_h0ij(1,1),tmp_h0ij(2,1)) 756 757 if (Mflags%has_vxc==1) & 758 & Mels%vxc(ib,jb,ik_ibz,is) = Mels%vxc(ib,jb,ik_ibz,is) + DCMPLX(tmp_xc(1,1),tmp_xc(2,1)) 759 760 if (Mflags%has_vxcval==1) & 761 & Mels%vxcval(ib,jb,ik_ibz,is) = Mels%vxcval(ib,jb,ik_ibz,is) + DCMPLX(tmp_xcval(1,1),tmp_xcval(2,1)) 762 763 if (Mflags%has_vxcval_hybrid==1) & 764 & Mels%vxcval_hybrid(ib,jb,ik_ibz,is) = Mels%vxcval_hybrid(ib,jb,ik_ibz,is) + DCMPLX(tmp_xcval(1,1),tmp_xcval(2,1)) 765 766 if (Mflags%has_vhartree==1) & 767 & Mels%vhartree(ib,jb,ik_ibz,is) = Mels%vhartree(ib,jb,ik_ibz,is) + DCMPLX(tmp_H (1,1),tmp_H (2,1)) 768 769 if (Mflags%has_vu==1) & 770 & Mels%vu(ib,jb,ik_ibz,is) = DCMPLX(tmp_U(1,1),tmp_U(2,1)) 771 772 if (Mflags%has_sxcore==1) & 773 & Mels%sxcore(ib,jb,ik_ibz,is) = DCMPLX(tmp_sigcx(1,1),tmp_sigcx(2,1)) 774 775 else 776 777 if (Mflags%has_hbare==1) & 778 & Mels%hbare(ib,jb,ik_ibz,:) = Mels%hbare(ib,jb,ik_ibz,:) + DCMPLX(tmp_h0ij(1,:),tmp_h0ij(2,:)) 779 780 if (Mflags%has_vxc==1) & 781 & Mels%vxc(ib,jb,ik_ibz,:) = Mels%vxc(ib,jb,ik_ibz,:) + DCMPLX(tmp_xc(1,:),tmp_xc(2,:)) 782 783 if (Mflags%has_vxcval==1) & 784 & Mels%vxcval(ib,jb,ik_ibz,:) = Mels%vxcval(ib,jb,ik_ibz,:) + DCMPLX(tmp_xcval(1,:),tmp_xcval(2,:)) 785 786 if (Mflags%has_vxcval_hybrid==1) & 787 & Mels%vxcval_hybrid(ib,jb,ik_ibz,:) = Mels%vxcval_hybrid(ib,jb,ik_ibz,:) + DCMPLX(tmp_xcval(1,:),tmp_xcval(2,:)) 788 789 if (Mflags%has_vhartree==1) & 790 & Mels%vhartree(ib,jb,ik_ibz,:) = Mels%vhartree(ib,jb,ik_ibz,:) + DCMPLX(tmp_H (1,:),tmp_H (2,:)) 791 792 if (Mflags%has_vu==1) & 793 & Mels%vu(ib,jb,ik_ibz,:) = DCMPLX(tmp_U(1,:),tmp_U(2,:)) 794 end if 795 796 end do !ib 797 end do !jb 798 799 end do !is 800 end do !ikc 801 802 ABI_FREE(dimlmn) 803 call pawcprj_free(Cprj_b1ks) 804 ABI_FREE(Cprj_b1ks) 805 call pawcprj_free(Cprj_b2ks) 806 ABI_FREE(Cprj_b2ks) 807 end if !PAW 808 809 ABI_FREE(bbp_ks_distrb) 810 811 ! Sum up contributions on each node 812 ! Set the corresponding has_* flags to 2. 813 call Mels%mpisum(wfd%comm) 814 815 ! Reconstruct lower triangle. 816 call Mels%herm() 817 818 ABI_FREE(kcalc2ibz) 819 call destroy_mpi_enreg(MPI_enreg_seq) 820 821 DBG_EXIT("COLL") 822 823 end subroutine calc_vhxc_me
ABINIT/m_vhxc_me [ Modules ]
NAME
m_vhxc_me
FUNCTION
Evaluate the matrix elements of $v_H$ and $v_{xc}$ and $v_U$
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_vhxc_me 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xcdata 28 use libxc_functionals 29 use m_dtset 30 use m_distribfft 31 32 use defs_datatypes,only : pseudopotential_type 33 use defs_abitypes, only : MPI_type 34 use m_pawang, only : pawang_type 35 use m_pawtab, only : pawtab_type 36 use m_paw_an, only : paw_an_type 37 use m_paw_ij, only : paw_ij_type 38 use m_pawfgrtab, only : pawfgrtab_type 39 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 40 use m_paw_denpot, only : paw_mknewh0 41 use m_hide_blas, only : xdotc 42 use m_wfd, only : wfdgw_t, wave_t 43 use m_crystal, only : crystal_t 44 use m_melemts, only : melflags_t, melements_t 45 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 46 use m_kg, only : mkkin 47 use m_rhotoxc, only : rhotoxc 48 49 implicit none 50 51 private