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