TABLE OF CONTENTS


ABINIT/calc_vhxc_me [ Functions ]

[ Top ] [ 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.

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

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

ABINIT/m_vhxc_me [ Modules ]

[ Top ] [ Modules ]

NAME

 m_vhxc_me

FUNCTION

  Evaluate the matrix elements of $v_H$ and $v_{xc}$ and $v_U$

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_vhxc_me
28 
29  implicit none
30 
31  private