Line data Source code
1 : !!****m* ABINIT/m_prep_calc_ucrpa
2 : !! NAME
3 : !! m_prep_calc_ucrpa
4 : !!
5 : !! FUNCTION
6 : !! Prepare data for the calculation of U with the CRPA method: oscillators strenghs and k-points.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! INPUTS
15 : !!
16 : !! OUTPUT
17 : !!
18 : !! SOURCE
19 :
20 : #if defined HAVE_CONFIG_H
21 : #include "config.h"
22 : #endif
23 :
24 :
25 : #include "abi_common.h"
26 :
27 : MODULE m_prep_calc_ucrpa
28 :
29 : use defs_basis
30 : use m_abicore
31 : use m_gwdefs!, only : czero_gw, cone_gw, j_gw, sigparams_t
32 : use m_xmpi
33 : use m_defs_ptgroups
34 : use m_errors
35 :
36 : use defs_datatypes, only : pseudopotential_type
37 : use m_time, only : timab
38 : use m_hide_blas, only : xdotc
39 : use m_geometry, only : normv
40 : use m_crystal, only : crystal_t
41 : use m_fft_mesh, only : rotate_FFT_mesh
42 : use m_bz_mesh, only : kmesh_t, findqg0
43 : use m_gsphere, only : gsphere_t
44 : use m_io_tools, only : flush_unit, open_file
45 : use m_vcoul, only : vcoul_t
46 : use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
47 : use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
48 : use m_pawang, only : pawang_type
49 : use m_pawtab, only : pawtab_type
50 : use m_pawfgrtab, only : pawfgrtab_type
51 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
52 : use m_paw_nhat, only : pawmknhat_psipsi
53 : use m_paw_sym, only : paw_symcprj
54 : use m_wfd, only : wfd_t, wave_t
55 : use m_oscillators, only : rho_tw_g
56 : use m_esymm, only : esymm_t
57 : use m_read_plowannier, only : read_plowannier
58 : use m_plowannier, only : plowannier_type,operwan_realspace_type
59 : use m_ebands, only : ebands_t
60 :
61 : implicit none
62 :
63 : private
64 :
65 : public :: prep_calc_ucrpa
66 : !!***
67 :
68 : contains
69 :
70 : !!****f* ABINIT/prep_calc_ucrpa
71 : !! NAME
72 : !! prep_calc_ucrpa
73 : !!
74 : !! FUNCTION
75 : !! Prepare data for the calculation of U with the CRPA method: oscillators strenghs and k-points.
76 : !!
77 : !! COPYRIGHT
78 : !! Copyright (C) 1999-2026 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf,TApplencourt,BAmadon)
79 : !! This file is distributed under the terms of the
80 : !! GNU General Public License, see ~abinit/COPYING
81 : !! or http://www.gnu.org/copyleft/gpl.txt .
82 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
83 : !!
84 : !! INPUTS
85 : !! sigmak_ibz=Index of the k-point in the IBZ.
86 : !! minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
87 : !! Gsph_x<gsphere_t>= Info on the G-sphere used for Sigma_x
88 : !! %nsym=number of symmetry operations
89 : !! %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
90 : !! operation and G is one of the ng vectors in reciprocal space
91 : !! %timrev=2 if time-reversal symmetry is used, 1 otherwise
92 : !! %gvec(3,ng)=integer coordinates of each plane wave in reciprocal space
93 : !! ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
94 : !! Kmesh <kmesh_t>
95 : !! %nbz=Number of points in the BZ
96 : !! %nibz=Number of points in IBZ
97 : !! %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
98 : !! %kbz(3,nbz)=k-point coordinates, full Brillouin zone
99 : !! %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
100 : !! %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
101 : !! %ktabp(nbz)= phase factor associated to tnons
102 : !! gwx_ngfft(18)=Information about 3D FFT for the oscillator strengths, see ~abinit/doc/variables/vargs.htm#ngfft
103 : !! gwx_nfftot=number of points of the FFT grid for GW wavefunctions
104 : !! Vcp <vcoul_t datatype> containing information on the cutoff technique
105 : !! %vc_sqrt(npwx,nqibz)= square-root of the coulombian potential for q-points in the IBZ
106 : !! Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
107 : !! Pawang <type(pawang_type)>=paw angular mesh and related data
108 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
109 : !! %usepaw=1 for PAW, 0 for NC pseudopotentials.
110 : !! Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
111 : !! %ibz=q points where $\tilde\epsilon^{-1}$ has been computed
112 : !! %bz(3,nqbz)=coordinates of all q-points in BZ
113 : !! Sigp <sigparams_t> (see the definition of this structured datatype)
114 : !! Cryst<crystal_t>=Info on unit cell and symmetries
115 : !! %natom=number of atoms in unit cell
116 : !! %ucvol=unit cell volume
117 : !! %nsym=number of symmetry operations
118 : !! %typat(natom)=type of each atom
119 : !! much slower but it requires less memory
120 : !! QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
121 : !! eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin
122 : !! occ(Sigp%nbnds,Kmesh%nibz,nsppol)=occupation numbers, for each k point in IBZ, each band and spin
123 : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
124 : !! allQP_sym(%nkibz,%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the
125 : !! little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
126 : !! prtvol=Flags governing verbosity level.
127 : !!
128 : !! OUTPUT
129 : !!
130 : !! NOTES
131 : !! 1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.
132 : !!
133 : !! 2) On the symmetrization of Sigma matrix elements
134 : !! If Sk = k+G0 then M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G} (k,q)
135 : !! If -Sk = k+G0 then M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)
136 : !!
137 : !! Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
138 : !! that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
139 : !! similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
140 : !! has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.
141 : !!
142 : !! 3) the unitary transformation relating wavefunctions
143 : !! at symmetric k-points should be taken into account during the symmetrization
144 : !! of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
145 : !! it is possible to make an invariant by just including all the degenerate states and
146 : !! averaging the final results over the degenerate subset. Here we divide the states
147 : !! where the QP energies are required into complexes. Note however that this approach is not
148 : !! based on group theory, and it might lead to spurious results in case of accidental degeneracies.
149 : !!
150 : !! SOURCE
151 :
152 0 : subroutine prep_calc_ucrpa(sigmak_ibz,ikcalc,itypatcor,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Gsph_x,Vcp,Kmesh,Qmesh,lpawu,&
153 0 : & M1_q_m,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,&
154 0 : & Psps,Wfd,Wfdf,allQP_sym,gwx_ngfft,ngfftf,&
155 0 : & prtvol,pawcross,plowan_compute,rhot1_q_m,wanbz,rhot1)
156 :
157 : #ifndef HAVE_CRPA_OPTIM
158 : #ifdef FC_INTEL
159 : #warning "optimization of m_prec_calc_ucrpa is deactivated on intel fortran"
160 : !DEC$ NOOPTIMIZE
161 : #endif
162 : #endif
163 :
164 : !Arguments ------------------------------------
165 : !scalars
166 : integer,intent(in) :: sigmak_ibz,ikcalc,itypatcor,prtvol,lpawu,minbnd,maxbnd,pawcross,plowan_compute
167 : type(crystal_t),intent(in) :: Cryst
168 : type(ebands_t),target,intent(in) :: QP_BSt
169 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
170 : type(vcoul_t),intent(in) :: Vcp
171 : type(gsphere_t),intent(in) :: Gsph_x
172 : ! type(littlegroup_t),intent(in) :: Ltg_k
173 : type(Pseudopotential_type),intent(in) :: Psps
174 : type(sigparams_t),target,intent(in) :: Sigp
175 : type(pawang_type),intent(in) :: Pawang
176 : class(wfd_t),target,intent(inout) :: Wfd,Wfdf
177 : !arrays
178 : complex(dp), intent(out) :: rhot1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,sigp%npwx,Qmesh%nibz)
179 : complex(dp), intent(out) :: M1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,sigp%npwx,Qmesh%nibz)
180 : integer,intent(in) :: gwx_ngfft(18),ngfftf(18)
181 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
182 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
183 : type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
184 : type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
185 : type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
186 : type(plowannier_type),intent(in) :: wanbz
187 : type(operwan_realspace_type),target,intent(inout) :: rhot1(Sigp%npwx,Qmesh%nibz)
188 :
189 : !Local variables ------------------------------
190 : !scalars
191 : integer,parameter :: use_pawnhat=0,ider0=0,ndat1=1
192 : integer :: bandinf,bandsup
193 : integer :: gwcalctyp,izero,ib_sum,ib,ib1,ib2,ig,ig_rot,ii,iik,itim_q,i2
194 : integer :: ik_bz,ik_ibz,isym_q,iq_bz,iq_ibz,spin,isym,itypatcor_read,jb,iat
195 : integer :: jik,jk_bz,jk_ibz,lcor,m1,m3,nspinor,nsppol,ifft
196 : integer :: ibsp,dimcprj_gw
197 : integer :: spad
198 : integer :: comm
199 : integer :: ispinor1,ispinor3,isym_kgw,isym_ki,gwx_mgfft,use_padfft,use_padfftf,gwx_fftalga,gwx_fftalgb
200 : integer :: gwx_nfftot,nfftf,mgfftf,ierr
201 : integer :: nhat12_grdim
202 : integer :: iatom1,iatom2,il1,il2,im1,im2,ispinor2,pos1,pos2,wan_jb,wan_ib_sum,pwx
203 : real(dp) :: fact_sp,theta_mu_minus_esum,tol_empty,norm,weight
204 : complex(dp) :: ctmp,scprod,ph_mkgwt,ph_mkt,eikr
205 : logical :: iscompatibleFFT,q_is_gamma
206 : character(len=500) :: msg
207 : type(wave_t),pointer :: wave_sum, wave_jb
208 : !arrays
209 : integer :: g0(3),spinor_padx(2,4)
210 0 : integer,pointer :: igfftxg0(:),igfftfxg0(:)
211 0 : integer,allocatable :: gwx_gfft(:,:),gwx_gbound(:,:),gboundf(:,:)
212 0 : integer,allocatable :: ktabr(:,:),irottb(:,:),ktabrf(:,:)
213 : real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),qbz(3),q0(3),tsec(2)
214 : real(dp) :: spinrot_kbz(4),spinrot_kgw(4)
215 0 : real(dp),pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
216 0 : real(dp),allocatable :: nhat12(:,:,:),grnhat12(:,:,:,:)
217 0 : complex(gwp),allocatable :: vc_sqrt_qbz(:)
218 0 : complex(gwp),allocatable :: rhotwg_ki(:,:)
219 0 : complex(gwp),allocatable :: wfr_bdgw(:,:),wfr_sum(:)
220 0 : complex(gwp),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
221 0 : complex(gwp),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:)
222 0 : complex(gwp),pointer :: cg_jb(:),cg_sum(:)
223 : complex(dp) :: ovlp(2)
224 0 : complex(dp),allocatable :: coeffW_BZ(:,:,:,:,:,:)
225 0 : complex(dp),pointer :: ptr_rhot(:,:,:,:,:)
226 0 : logical :: can_symmetrize(Wfd%nsppol)
227 0 : logical,allocatable :: bks_mask(:,:,:)
228 0 : type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
229 0 : type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
230 0 : type(esymm_t),pointer :: QP_sym(:)
231 : !type(plowannier_type) :: wan
232 : logical :: ecriture=.FALSE.
233 : logical :: l_ucrpa,luwindow
234 : integer :: g0_dump(3),iq_ibz_dump,dumint(2)
235 :
236 : !************************************************************************
237 :
238 0 : l_ucrpa=.true.
239 :
240 : DBG_ENTER("COLL")
241 :
242 : !
243 : ! === Initial check ===
244 0 : ABI_CHECK(Sigp%npwx==Gsph_x%ng,'')
245 :
246 0 : call timab(430,1,tsec) ! csigme (SigX)
247 :
248 0 : gwcalctyp=Sigp%gwcalctyp
249 : !
250 : ! === Initialize MPI variables ===
251 0 : comm = Wfd%comm
252 :
253 : !
254 : ! === Initialize some values ===
255 0 : nspinor = Wfd%nspinor
256 0 : nsppol = Wfd%nsppol
257 0 : spinor_padx(:,:)=RESHAPE((/0,0,Sigp%npwx,Sigp%npwx,0,Sigp%npwx,Sigp%npwx,0/),(/2,4/))
258 :
259 0 : qp_ene => QP_BSt%eig(:,:,:)
260 0 : qp_occ => QP_BSt%occ(:,:,:)
261 :
262 : ! Exctract the symmetries of the bands for this k-point
263 0 : QP_sym => allQP_sym(sigmak_ibz,1:nsppol)
264 :
265 0 : ib1=minbnd
266 0 : ib2=maxbnd
267 :
268 : ! === Read Wannier function coefficients for Ucrpa
269 : ! === for future computation of rhot_m_q directly in this routine.
270 :
271 0 : dumint=0
272 0 : luwindow=.true.
273 : ! write(6,*) "cc",allocated(coeffW_BZ)
274 0 : if (plowan_compute <10)then
275 : call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor_read,Kmesh,lcor,luwindow,&
276 0 : & nspinor,nsppol,pawang,prtvol,dumint)
277 0 : if(lcor/=lpawu) then
278 0 : msg = "lcor and lpawu differ in prep_calc_ucrpa"
279 0 : ABI_ERROR(msg)
280 : endif
281 : endif
282 :
283 : ! === End of read Wannier function coefficients for Ucrpa
284 :
285 :
286 : !
287 : ! === Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
288 0 : jk_bz=Sigp%kptgw2bz(ikcalc)
289 : !write(6,*) "ikcalc,jk_bz",ikcalc,jk_bz
290 : !write(6,*) "ikcalc",Kmesh%bz(:,ikcalc)
291 : !write(6,*) "jk_bz",Kmesh%bz(:,jk_bz)
292 : ! jk_bz=ikcalc
293 0 : call kmesh%get_BZ_item(jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
294 : ! write(6,*) "jk_ibz",Kmesh%ibz(:,jk_ibz)
295 : ! write(6,*) "jk_bz,jk_ibz",jk_bz,jk_ibz,isym_kgw,itim
296 : !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
297 0 : spinrot_kgw(:)=Cryst%spinrot(:,isym_kgw)
298 : !
299 0 : write(msg,'(2a,3f8.3,a,i4,a,2(i3,a))')ch10,&
300 0 : & ' Calculating Oscillator element at k= ',kgw, "k-point number",ikcalc,&
301 0 : & ' bands n = from ',ib1,' to ',ib2,ch10
302 0 : call wrtout(std_out,msg,'COLL')
303 :
304 :
305 0 : if (ANY(gwx_ngfft(1:3) /= Wfd%ngfft(1:3)) ) then
306 0 : call wfd%change_ngfft(Cryst,Psps,gwx_ngfft)
307 : end if
308 0 : gwx_mgfft = MAXVAL(gwx_ngfft(1:3))
309 0 : gwx_fftalga = gwx_ngfft(7)/100
310 0 : gwx_fftalgb = MOD(gwx_ngfft(7),100)/10
311 :
312 0 : if (pawcross==1) then
313 0 : mgfftf = MAXVAL(ngfftf(1:3))
314 : end if
315 :
316 0 : can_symmetrize = .FALSE.
317 0 : if (Sigp%symsigma>0) then
318 0 : can_symmetrize = .TRUE.
319 0 : if (gwcalctyp >= 20) then
320 0 : do spin=1,Wfd%nsppol
321 0 : can_symmetrize(spin) = .not. QP_sym(spin)%failed()
322 0 : if (.not.can_symmetrize(spin)) then
323 : write(msg,'(a,i0,4a)')&
324 0 : " Symmetrization cannot be performed for spin: ",spin,ch10,&
325 0 : " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
326 0 : ABI_WARNING(msg)
327 : end if
328 : end do
329 : end if
330 0 : ABI_CHECK(nspinor==1,'Symmetrization with nspinor=2 not implemented')
331 : end if
332 :
333 0 : ABI_MALLOC(rhotwg_ki,(Sigp%npwx*nspinor,minbnd:maxbnd))
334 0 : rhotwg_ki=czero_gw
335 0 : ABI_MALLOC(vc_sqrt_qbz,(Sigp%npwx))
336 : !
337 : ! === Normalization of theta_mu_minus_esum ===
338 : ! * If nsppol==2, qp_occ $\in [0,1]$
339 0 : SELECT CASE (nsppol)
340 : CASE (1)
341 0 : fact_sp=half; tol_empty=0.01 ! below this value the state is assumed empty
342 0 : if (Sigp%nspinor==2) then
343 0 : fact_sp=one; tol_empty=0.005 ! below this value the state is assumed empty
344 : end if
345 : CASE (2)
346 0 : fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic
347 : CASE DEFAULT ! spin unpolarized system is treated using nsppol==2
348 0 : ABI_BUG('Wrong nsppol')
349 : END SELECT
350 :
351 : ! Remove empty states from the list of states that will be distributed.
352 0 : ABI_MALLOC(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol))
353 0 : bks_mask=.FALSE.
354 0 : do spin=1,nsppol
355 0 : do ik_bz=1,Kmesh%nbz
356 0 : ik_ibz = Kmesh%tab(ik_bz)
357 0 : do ib_sum=1,Sigp%nbnds
358 0 : bks_mask(ib_sum,ik_bz,spin) = (qp_occ(ib_sum,ik_ibz,spin)>=tol_empty)
359 : end do
360 : end do
361 : end do
362 :
363 : ! ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol))
364 : ! call sigma_distribution(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
365 : ! call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
366 : ! write(6,*)"lim", ib1,ib2
367 : ! do ib_sum= ib1,ib2
368 : ! do ik_bz=1, Kmesh%nbz
369 : ! write(6,*) ib_sum,ik_bz, proc_distrb(ib_sum,ik_bz,1),Wfd%my_rank
370 : ! enddo
371 : ! enddo
372 :
373 0 : ABI_FREE(bks_mask)
374 :
375 0 : write(msg,'(a,i8)')" Will sum all (b,k,s) occupied states in Sigma_x for k-point",ikcalc
376 0 : call wrtout(std_out,msg,'PERS')
377 : !
378 : ! The index of G-G0 in the FFT mesh the oscillators ===
379 : ! * Sigp%mG0 gives the MAX G0 component to account for umklapp.
380 : ! * Note the size MAX(Sigp%npwx,Sigp%npwc).
381 0 : ABI_MALLOC(igfftxg0,(Gsph_x%ng))
382 : !
383 : ! === Precalculate the FFT index of $ R^{-1}(r-\tau) $ ===
384 : ! * S=\transpose R^{-1} and k_BZ = S k_IBZ
385 : ! * irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
386 0 : gwx_nfftot = PRODUCT(gwx_ngfft(1:3))
387 0 : ABI_MALLOC(irottb,(gwx_nfftot,Cryst%nsym))
388 0 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwx_ngfft,irottb,iscompatibleFFT)
389 0 : if (.not.iscompatibleFFT) then
390 0 : msg = "FFT mesh is not compatible with symmetries. Results might be affected by large errors!"
391 0 : ABI_WARNING(msg)
392 : end if
393 :
394 0 : ABI_MALLOC(ktabr,(gwx_nfftot,Kmesh%nbz))
395 0 : do ik_bz=1,Kmesh%nbz
396 0 : isym=Kmesh%tabo(ik_bz)
397 0 : do ifft=1,gwx_nfftot
398 0 : ktabr(ifft,ik_bz)=irottb(ifft,isym)
399 : end do
400 : end do
401 0 : ABI_FREE(irottb)
402 :
403 0 : if (Psps%usepaw==1 .and. pawcross==1) then
404 0 : nfftf = PRODUCT(ngfftf(1:3))
405 0 : ABI_MALLOC(irottb,(nfftf,Cryst%nsym))
406 0 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT)
407 :
408 0 : ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz))
409 0 : do ik_bz=1,Kmesh%nbz
410 0 : isym=Kmesh%tabo(ik_bz)
411 0 : do ifft=1,nfftf
412 0 : ktabrf(ifft,ik_bz)=irottb(ifft,isym)
413 : end do
414 : end do
415 0 : ABI_FREE(irottb)
416 : end if
417 : !
418 : ! === Additional allocations for PAW ===
419 0 : if (Psps%usepaw==1) then
420 0 : ABI_MALLOC(Cprj_ksum,(Cryst%natom,nspinor))
421 0 : call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
422 :
423 0 : nhat12_grdim=0
424 : if (use_pawnhat==1) then ! Compensation charge for \phi_a^*\phi_b
425 : call wrtout(std_out,"Using nhat12","COLL")
426 : ABI_MALLOC(nhat12 ,(2,gwx_nfftot,nspinor**2))
427 : ABI_MALLOC(grnhat12,(2,gwx_nfftot,nspinor**2,3*nhat12_grdim))
428 : end if
429 : end if ! usepaw==1
430 : !
431 :
432 : if (Sigp%symsigma>0) then
433 : !call littlegroup_print(Ltg_k,std_out,prtvol,'COLL')
434 : !
435 : ! === Find number of complexes and number of bands in each complex ===
436 : ! The tolerance is a little bit arbitrary (0.001 eV)
437 : ! It could be reduced, in particular in case of nearly accidental degeneracies
438 : ! if (ANY(degtab/=0)) then ! If two states do not belong to the same complex => matrix elements of v_xc differ
439 : ! write(msg,'(a,3f8.3,a)')' Degenerate states at k-point = ( ',kgw(:),' ).'
440 : ! call wrtout(std_out,msg,'COLL')
441 : ! do spin=1,nsppol
442 : ! do ib=ib1,ib2
443 : ! do jb=ib+1,ib2
444 : ! if (degtab(ib,jb,spin)==1) then
445 : ! write(msg,'(a,i2,a,i4,a,i4)')' (spin ',spin,')',ib,' <====> ',jb
446 : ! call wrtout(std_out,msg,'COLL')
447 : ! if (ABS(Sr%vxcme(ib,jk_ibz,spin)-Sr%vxcme(jb,jk_ibz,spin))>ABS(tol6*Sr%vxcme(jb,jk_ibz,spin))) then
448 : ! write(msg,'(7a)')&
449 : !& ' It seems that an accidental degeneracy is occurring at this k-point ',ch10,&
450 : !& ' In this case, using symsigma=1 might lead to spurious results as the algorithm ',ch10,&
451 : !& ' will treat these states as degenerate, and it won''t be able to remove the degeneracy. ',ch10,&
452 : !& ' In order to avoid this deficiency, run the calculation using symsigma=0'
453 : ! ABI_WARNING(msg)
454 : ! end if
455 : ! end if
456 : ! end do
457 : ! end do
458 : ! end do
459 : ! end if
460 : end if !symsigma
461 :
462 0 : ABI_MALLOC(wfr_sum,(gwx_nfftot*nspinor))
463 0 : if (pawcross==1) then
464 0 : ABI_MALLOC(ur_ae_sum,(nfftf*nspinor))
465 0 : ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor))
466 0 : ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor))
467 : end if
468 :
469 : !!*******************************************
470 : !! Save ik_bz and Norm of G vectors.
471 : !! FOR THE UCRPA calculation
472 : !!*******************************************
473 0 : if (ikcalc==1) then
474 0 : ecriture=.TRUE.
475 : ! open(unit=2011,file='ikbz_COORD',form='formatted',status='unknown')
476 : ! do ik_bz=1,Kmesh%nbz
477 : ! call get_BZ_item(Kmesh,ik_bz,&
478 : ! kbz_coord,ik_ibz,isym_kgw,iik,ph_mkt)
479 : ! write(2011,*) ik_bz,kbz_coord(:)
480 : ! end do
481 : ! close(2011)
482 :
483 : ! if (prtvol>10.and.jk_bz==1) then ! probably just to print one time.
484 : ! !!q=0 Forcement donc divergence pour G=0
485 : ! if (open_file("normeG", msg, unit=2022, form="formatted", status="unknown") /= 0) then
486 : ! ABI_ERROR(msg)
487 : ! end if
488 : ! write(2022,*) 1,real(CMPLX(SQRT(Vcp%i_sz),0.0_gwp)),real((4*3.14159265)**(0.5)/CMPLX(SQRT(Vcp%i_sz),0.0_gwp))
489 : ! write G=0 term for the potential computed elsewhere.
490 : !
491 : ! do ig=2,Sigp%npwx
492 : ! write(2022,*) ig,real(Vcp%vc_sqrt(ig,1)),real((4*3.14159265)**(0.5)/Vcp%vc_sqrt(ig,1))
493 : ! end do
494 : ! write potential for ig, q=0
495 : ! close(2022)
496 : ! end if
497 :
498 : ! write header for q point written later.
499 : ! if (Wfd%my_rank==0) then
500 : ! open(unit=2015,file='iqbz_COORD',form='formatted',status='unknown')
501 : ! write (2015,*) "q pour le k", jk_bz, ikcalc
502 : ! close(2015)
503 : ! endif
504 :
505 : else
506 0 : ecriture=.FALSE.
507 : end if
508 :
509 : !!*******************************************
510 : !! End of print if ik_bz for UCRPA calc
511 : !!*******************************************
512 :
513 : !
514 : ! =======================================
515 : ! ==== Begin loop over k_i in the BZ ====
516 : ! =======================================
517 :
518 0 : do spin=1,nsppol
519 :
520 : ! if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE
521 : ! write(6,*) "AA",Wfd%my_rank,spin
522 : !
523 : ! * Load wavefunctions for GW corrections.
524 0 : ABI_STAT_MALLOC(wfr_bdgw,(gwx_nfftot*nspinor,ib1:ib2), ierr)
525 0 : ABI_CHECK(ierr==0, "out of memory in wfr_bdgw")
526 0 : do jb=ib1,ib2
527 0 : call wfd%get_ur(jb,jk_ibz,spin,wfr_bdgw(:,jb))
528 : ! write(6,'(a,6i4)')"indforwfd" ,jb,jk_ibz,spin
529 : end do
530 :
531 0 : if (Wfd%usepaw==1) then ! * Load cprj for GW states, note the indexing.
532 0 : dimcprj_gw=nspinor*(ib2-ib1+1)
533 0 : ABI_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
534 0 : call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
535 0 : ibsp=ib1
536 0 : do jb=ib1,ib2
537 0 : call wfd%get_cprj(jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
538 0 : call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
539 0 : call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
540 0 : ibsp=ibsp+nspinor
541 : end do
542 0 : if (pawcross==1) then
543 0 : ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,ib1:ib2))
544 0 : ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,ib1:ib2))
545 0 : ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,ib1:ib2))
546 0 : do jb=ib1,ib2
547 : call wfdf%paw_get_aeur(jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
548 0 : & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
549 0 : ur_ae_bdgw(:,jb)=ur_ae_sum
550 0 : ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum
551 0 : ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum
552 : end do
553 : end if
554 : end if
555 :
556 0 : do ik_bz=1,Kmesh%nbz
557 : !
558 : ! === Parallelization over k-points and spin ===
559 : ! * For the spin there is another check in the inner loop
560 : ! if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
561 : ! write(6,*) "BB",Wfd%my_rank,spin,ik_bz
562 : !
563 : ! * Find the corresponding irreducible k-point
564 0 : call kmesh%get_BZ_item(ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt)
565 0 : spinrot_kbz(:)=Cryst%spinrot(:,isym_ki)
566 : ! write(6,'(a,6i4)')"indices" ,jk_bz,jk_ibz,ik_bz,ik_ibz,spin
567 :
568 : ! * Identify q and G0 where q+G0=k_GW-k_i
569 0 : kgw_m_ksum=kgw-ksum
570 : ! write(6,*) "kgw ",kgw
571 : ! write(6,*) "ksum ",ksum
572 0 : call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0)
573 :
574 : ! if(iq_bz/=1.or.ik_bz/=1) cycle
575 : ! write(6,*) "g0",g0
576 : ! write(6,*) " ik_bz=",ik_bz
577 : ! write(6,*) " iq_bz=",iq_bz
578 :
579 : ! === Symmetrize the matrix elements ===
580 : ! * Sum only q"s in IBZ_k. In this case elements are weighted
581 : ! according to wtqp and wtqm. wtqm is for time-reversal.
582 : ! wtqp=1; wtqm=0
583 : ! if (can_symmetrize(spin)) then
584 : ! ! if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
585 : ! wtqp=0; wtqm=0
586 : ! do isym=1,Ltg_k%nsym_sg
587 : ! wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz)
588 : ! wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz)
589 : ! end do
590 : ! end if
591 :
592 :
593 : !
594 : ! * Find the corresponding irreducible q-point.
595 0 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
596 0 : q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOLQ0)
597 :
598 : !!*******************************************
599 : !! Check if qbz belongs to IBZ.
600 : !! because dielectric matrix is computed in the IBZ.
601 : !! then output the q point.
602 : !!*******************************************
603 : ! write(6,*) "kkk1",ik_bz,jk_bz,iq_ibz
604 :
605 0 : if (.NOT. qmesh%has_IBZ_item(qbz,iq_ibz_dump,g0_dump)) then
606 : cycle
607 : end if
608 :
609 0 : write(msg,'(2(a,i4),a,i3)')' prep_calc_ucrpa : ik_bz ',ik_bz,'/',Kmesh%nbz,' done'
610 0 : call wrtout(std_out,msg,'PERS')
611 : ! write(6,*) "kkk1p",ik_bz,jk_bz,iq_ibz
612 :
613 : ! write(std_out,*)'prep_calc_ucrpa:ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank,"iq_BZ",iq_bz,"iQ_iBZ",iq_ibz
614 : ! write(123,*)'prep_calc_ucrpa:jk_bz ',jk_bz,'ikmq_bz',ik_bz,'iq_ibz',iq_ibz
615 :
616 : !Ecriture du iq_ibz
617 : ! if (ecriture.and.Wfd%my_rank==0) then
618 : ! open(unit=2016,file='iqbz_COORD',form='formatted',status='unknown',position='append')
619 : ! write(2016,*) iq_ibz,qbz(:),Qmesh%wt(iq_ibz)
620 : ! ! write(2011,*) iq_bz,qbz(:)
621 : ! close(2016)
622 : ! end if
623 : !!*******************************************
624 : !! End of modif for UCRPA.
625 : !!*******************************************
626 :
627 : !
628 : ! Tables for the FFT of the oscillators.
629 : ! a) FFT index of the G-G0.
630 0 : ABI_MALLOC(gwx_gbound,(2*gwx_mgfft+8,2))
631 0 : call Gsph_x%fft_tabs(g0,gwx_mgfft,gwx_ngfft,use_padfft,gwx_gbound,igfftxg0)
632 :
633 0 : if ( ANY(gwx_fftalga == (/2,4/)) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
634 : #ifdef FC_IBM
635 : ! XLF does not deserve this optimization (problem with [v67mbpt][t03])
636 : use_padfft = 0
637 : #endif
638 0 : if (use_padfft==0) then
639 0 : ABI_FREE(gwx_gbound)
640 0 : ABI_MALLOC(gwx_gbound,(2*gwx_mgfft+8,2*use_padfft))
641 : end if
642 :
643 0 : if (pawcross==1) then
644 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2))
645 0 : ABI_MALLOC(igfftfxg0,(Gsph_x%ng))
646 0 : call Gsph_x%fft_tabs(g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftfxg0)
647 0 : if ( ANY(gwx_fftalga == (/2,4/)) ) use_padfftf=0
648 0 : if (use_padfftf==0) then
649 0 : ABI_FREE(gboundf)
650 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
651 : end if
652 : end if
653 : !
654 : ! === Evaluate oscillator matrix elements ===
655 : ! * $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form.
656 0 : if (Psps%usepaw==1.and.use_pawnhat==0) then
657 0 : q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
658 0 : ABI_MALLOC(Pwij_qg,(Psps%ntypat))
659 0 : call pawpwij_init(Pwij_qg,Sigp%npwx,q0,Gsph_x%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
660 : end if
661 : !
662 : ! === Get Fourier components of the Coulomb interaction in the BZ ===
663 : ! * In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G|
664 : ! * The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
665 0 : do ig=1,Sigp%npwx
666 0 : ig_rot = Gsph_x%rottb(ig,itim_q,isym_q)
667 0 : vc_sqrt_qbz(ig_rot)=Vcp%vc_sqrt(ig,iq_ibz)
668 : end do
669 : !
670 : ! write(6,*) "kkk2",ik_bz,jk_bz,iq_ibz
671 : ! === Sum over bands ===
672 : !do ib_sum=1,Sigp%nbnds
673 0 : do ib_sum=ib1,ib2
674 : !write(6,*) "ib_sum",ib_sum
675 : !
676 : ! === Parallelism over spin ===
677 : ! * This processor has this k-point but what about spin?
678 : ! if (proc_distrb(ib_sum,ik_bz,spin)/=Wfd%my_rank) CYCLE
679 : ! write(6,*) "CC",Wfd%my_rank,spin,ik_bz,ib_sum
680 : !
681 : ! * Skip empty states.
682 : !if (qp_occ(ib_sum,ik_ibz,spin)<tol_empty) CYCLE
683 :
684 0 : call wfd%get_ur(ib_sum,ik_ibz,spin,wfr_sum)
685 : ! write(6,'(a,3i4)')"indforwfd2" ,ib_sum,ik_ibz,spin
686 : ! write(6,*) wfd_ihave_ug(Wfd,ib_sum,ik_ibz,spin,"Stored")
687 : ! write(6,*) wfd_ihave_ur(Wfd,ib_sum,ik_ibz,spin,"Stored")
688 :
689 0 : if (Psps%usepaw==1) then ! Load cprj for point ksum, this spin or spinor and *THIS* band.
690 : ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
691 : ! For sure there is a better way to do this!
692 0 : call wfd%get_cprj(ib_sum,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
693 0 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
694 0 : if (pawcross==1) then
695 : call wfdf%paw_get_aeur(ib_sum,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
696 0 : & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
697 : end if
698 : end if
699 :
700 0 : do jb=ib1,ib2 ! Get all <k-q,ib_sum,s|e^{-i(q+G).r}|s,jb,k>
701 : ! if(ib_sum.ne.jb) cycle
702 : ! if(ib_sum.ne.1.and.ib_sum.ne.10) cycle
703 : ! write(6,*) "jb",jb
704 :
705 : if (Psps%usepaw==1.and.use_pawnhat==1) then
706 : i2=jb; if (nspinor==2) i2=(2*jb-1)
707 : spad=(nspinor-1)
708 :
709 : izero=0
710 : call pawmknhat_psipsi(Cprj_ksum,Cprj_kgw(:,i2:i2+spad),ider0,izero,Cryst%natom,&
711 : & Cryst%natom,gwx_nfftot,gwx_ngfft,nhat12_grdim,nspinor,Cryst%ntypat,1,1,Pawang,Pawfgrtab,&
712 : & grnhat12,nhat12,pawtab)
713 :
714 : #if 1
715 : msg = "reinstate optional Argument in rho_tw_g but mind inca slave!"
716 : ABI_ERROR(msg)
717 : #else
718 : call rho_tw_g(nspinor,Sigp%npwx,gwx_nfftot,ndat1,gwx_ngfft,1,use_padfft,igfftxg0,gwx_gbound,&
719 : & wfr_sum ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz,&
720 : & wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
721 : & nspinor,rhotwg_ki(:,jb),nhat12=nhat12)
722 : #endif
723 :
724 : else
725 : call rho_tw_g(nspinor,Sigp%npwx,gwx_nfftot,ndat1,gwx_ngfft,1,use_padfft,igfftxg0,gwx_gbound,&
726 : & wfr_sum ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz,&
727 : & wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
728 0 : & nspinor,rhotwg_ki(:,jb))
729 :
730 0 : if (Psps%usepaw==1.and.use_pawnhat==0) then ! Add on-site contribution, projectors are already in BZ.
731 0 : i2=jb; if (nspinor==2) i2=(2*jb-1)
732 0 : spad=(nspinor-1)
733 : call paw_rho_tw_g(cryst, Pwij_qg,Sigp%npwx,nspinor,nspinor,Gsph_x%gvec,&
734 0 : Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb))
735 : end if
736 0 : if(iq_bz==1) then
737 0 : if((ib_sum/=jb).and.(abs(rhotwg_ki(1,jb))>tol8)) then
738 0 : if((ib_sum/=jb).and.(abs(rhotwg_ki(1,jb))>0.01_dp)) then
739 0 : write(std_out,*) "Warning: precision is low, oscillator strengh should be zero and is :",rhotwg_ki(1,jb)
740 : !else
741 : ! write(std_out,*) "Warning1: oscillator strengh",rhotwg_ki(1,jb)
742 : endif
743 : endif
744 0 : if((ib_sum==jb).and.(abs(rhotwg_ki(1,jb)-1_dp)>tol8)) then
745 0 : if((ib_sum==jb).and.(abs(rhotwg_ki(1,jb)-1_dp)>0.01_dp)) then
746 0 : write(std_out,*) "Warning: precision is low, oscillator strengh should be one and is :",rhotwg_ki(1,jb)
747 : !else
748 : ! write(std_out,*) "Warning1: oscillator strengh",rhotwg_ki(1,jb)
749 : endif
750 : endif
751 : endif
752 0 : if (Psps%usepaw==1.and.pawcross==1) then ! Add paw cross term
753 : call paw_cross_rho_tw_g(nspinor,Sigp%npwx,nfftf,ngfftf,1,use_padfftf,igfftfxg0,gboundf,&
754 : & ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,&
755 : & ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
756 0 : & nspinor,rhotwg_ki(:,jb))
757 : end if
758 : end if
759 :
760 : ! ************************************8
761 : ! Compute M Matrix in Wannier basis
762 : ! ************************************8
763 :
764 0 : if (ib_sum.GE.ib1.AND.ib_sum.LE.ib2) then
765 0 : call flush_unit(std_out)
766 0 : call flush_unit(ab_out)
767 0 : if (plowan_compute <10)then
768 0 : do iat=1, cryst%nattyp(itypatcor)
769 0 : do ispinor1=1,nspinor
770 0 : do ispinor3=1,nspinor
771 0 : do m1=1,2*lcor+1
772 0 : do m3=1,2*lcor+1
773 : M1_q_m(iat,ispinor1,ispinor3,m1,m3,:,iq_ibz)=M1_q_m(iat,ispinor1,ispinor3,m1,m3,:,iq_ibz)+&
774 0 : & rhotwg_ki(:,jb)*coeffW_BZ(iat,spin,jb,jk_bz,ispinor3,m3)*conjg(coeffW_BZ(iat,spin,ib_sum,ik_bz,ispinor1,m1))
775 : enddo
776 : enddo
777 : enddo
778 : enddo
779 : enddo
780 : end if
781 : endif
782 : ! ************************************8
783 : ! ************************************8
784 : !
785 : ! === Multiply by the square root of the Coulomb term ===
786 : ! * In 3-D systems, the factor sqrt(4pi) is included)
787 0 : do ii=1,nspinor
788 0 : spad=(ii-1)*Sigp%npwx
789 : !!$omp parallel workshare
790 0 : rhotwg_ki(spad+1:spad+Sigp%npwx,jb)=rhotwg_ki(spad+1:spad+Sigp%npwx,jb)*vc_sqrt_qbz(1:Sigp%npwx)
791 : !!$omp end parallel workshare
792 : end do
793 : !
794 : ! === Treat analytically the case q --> 0 ===
795 : ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
796 : ! while the Colulomb term is integrated out.
797 : ! * In the scalar case we have nonzero contribution only if ib_sum==jb
798 : ! * For nspinor==2 evalute <ib_sum,up|jb,up> and <ib_sum,dwn|jb,dwn>,
799 : ! impose orthonormalization since npwwfn might be < npwvec.
800 0 : if (ik_bz==jk_bz) then
801 0 : if (nspinor==1) then
802 0 : rhotwg_ki(1,jb)=czero_gw
803 0 : if (ib_sum==jb) rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)
804 : else
805 : ! TODO Recheck this!
806 0 : ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave_sum, msg) == 0, msg)
807 0 : cg_sum => wave_sum%ug
808 0 : ABI_CHECK(wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg)
809 0 : cg_jb => wave_jb%ug
810 :
811 0 : ctmp = xdotc(Wfd%npwarr(ik_ibz)*Wfd%nspinor,cg_sum,1,cg_jb,1)
812 0 : ovlp(1) = REAL(ctmp)
813 0 : ovlp(2) = AIMAG(ctmp)
814 :
815 0 : if (Psps%usepaw==1) then
816 0 : i2=(2*jb-1)
817 0 : ovlp = ovlp + paw_overlap(Cprj_ksum,Cprj_kgw(:,i2:i2+1),Cryst%typat,Pawtab)
818 : end if
819 : !ovlp(2) = -ovlp(1)
820 : !if (ib_sum==jb) ovlp(2)=cone_gw-ovlp(1)
821 0 : if (ib_sum==jb) then
822 0 : norm=DBLE(ovlp(1)+ovlp(2))
823 0 : ovlp(1)=DBLE(ovlp(1)/norm)
824 0 : ovlp(2)=DBLE(ovlp(2)/norm)
825 : else
826 0 : scprod=ovlp(1)+ovlp(2)
827 0 : ovlp(1)=ovlp(1)-scprod*half
828 0 : ovlp(2)=ovlp(2)-scprod*half
829 : end if
830 0 : rhotwg_ki(1 ,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)*ovlp(1)
831 0 : rhotwg_ki(Sigp%npwx+1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)*ovlp(2)
832 : end if
833 : end if
834 :
835 : ! ************************************8
836 : ! compute rhotwidle in Wannier basis for UcRPA (summed over k)
837 : ! ************************************8
838 : ! if (cryst%nsym==1) then
839 0 : weight=one
840 : ! else
841 : ! weight=Kmesh%wt(jk_bz)
842 : ! endif
843 0 : if (ib_sum.GE.ib1.AND.ib_sum.LE.ib2) then
844 0 : call flush_unit(std_out)
845 0 : call flush_unit(ab_out)
846 0 : if (plowan_compute<10)then
847 0 : do iat=1, cryst%nattyp(itypatcor)
848 0 : do ispinor1=1,nspinor
849 0 : do ispinor3=1,nspinor
850 0 : do m1=1,2*lcor+1
851 0 : do m3=1,2*lcor+1
852 : if(m1==2.and.m3==2) then
853 : endif
854 : rhot1_q_m(iat,ispinor1,ispinor3,m1,m3,:,iq_ibz)=&
855 : &rhot1_q_m(iat,ispinor1,ispinor3,m1,m3,:,iq_ibz)+&
856 : &rhotwg_ki(:,jb)*coeffW_BZ(iat,spin,jb,jk_bz,ispinor3,m3)&
857 0 : &*conjg(coeffW_BZ(iat,spin,ib_sum,ik_bz,ispinor1,m1))*weight
858 : enddo
859 : enddo
860 : enddo
861 : enddo
862 : enddo
863 : else
864 0 : wan_jb=jb-wanbz%bandi_wan+1
865 0 : wan_ib_sum=ib_sum-wanbz%bandi_wan+1
866 0 : do pwx=1,sigp%npwx
867 0 : do iatom1=1,wanbz%natom_wan
868 0 : do iatom2=1,wanbz%natom_wan
869 : !Loig Vaugier PhD eq. 5.11
870 : eikr=exp(- cmplx(0.0,1.0) * two_pi * ( &
871 : kmesh%bz(1,ik_bz)* ( cryst%xred(1,wanbz%iatom_wan(iatom1)) - cryst%xred(1,wanbz%iatom_wan(iatom2)) )+&
872 : kmesh%bz(2,ik_bz)* ( cryst%xred(2,wanbz%iatom_wan(iatom1)) - cryst%xred(2,wanbz%iatom_wan(iatom2)) )+&
873 0 : kmesh%bz(3,ik_bz)* ( cryst%xred(3,wanbz%iatom_wan(iatom1)) - cryst%xred(3,wanbz%iatom_wan(iatom2)) )))
874 0 : do pos1=1,size(wanbz%nposition(iatom1)%pos,1)
875 0 : do pos2=1,size(wanbz%nposition(iatom2)%pos,1)
876 0 : do il1=1,wanbz%nbl_atom_wan(iatom1)
877 0 : do il2=1,wanbz%nbl_atom_wan(iatom2)
878 0 : ptr_rhot=>rhot1(pwx,iq_ibz)%atom_index(iatom1,iatom2)%position(pos1,pos2)%atom(il1,il2)%matl
879 0 : do im1=1,2*wanbz%latom_wan(iatom1)%lcalc(il1)+1
880 0 : do im2=1,2*wanbz%latom_wan(iatom2)%lcalc(il2)+1
881 0 : do ispinor1=1,wanbz%nspinor
882 0 : do ispinor2=1,wanbz%nspinor
883 : ptr_rhot(im1,im2,spin,ispinor1,ispinor2)=&
884 : &ptr_rhot(im1,im2,spin,ispinor1,ispinor2)+&
885 : &rhotwg_ki(pwx,jb)*wanbz%psichi(jk_bz,wan_jb,iatom1)%atom(il1)%matl(im1,spin,ispinor1)*&
886 : &conjg(wanbz%psichi(ik_bz,wan_ib_sum,iatom2)%atom(il2)%matl(im2,spin,ispinor2))*weight&
887 0 : *eikr
888 : enddo!im2
889 : enddo!im1
890 : enddo!il2
891 : enddo!il1
892 : enddo!pos2
893 : enddo!pos1
894 : enddo!iatom2
895 : enddo!iatom1
896 : enddo!ispinor2
897 : enddo!ispinor1
898 : enddo!pwx
899 : endif!plowan_compute<10
900 : end if
901 : ! ************************************8
902 : ! ************************************8
903 :
904 : end do !jb Got all matrix elements from minbnd up to maxbnd.
905 :
906 0 : theta_mu_minus_esum=fact_sp*qp_occ(ib_sum,ik_ibz,spin)
907 :
908 : end do !ib_sum
909 : !
910 : ! Deallocate k-dependent quantities.
911 0 : ABI_FREE(gwx_gbound)
912 0 : if (pawcross==1) then
913 0 : ABI_FREE(gboundf)
914 : end if
915 :
916 0 : if (Psps%usepaw==1.and.use_pawnhat==0) then
917 0 : call pawpwij_free(Pwij_qg)
918 0 : ABI_FREE(Pwij_qg)
919 : end if
920 :
921 : end do !ik_bz Got all diagonal (off-diagonal) matrix elements.
922 :
923 0 : ABI_FREE(wfr_bdgw)
924 0 : if (Wfd%usepaw==1) then
925 0 : call pawcprj_free(Cprj_kgw )
926 0 : ABI_FREE(Cprj_kgw)
927 0 : if (pawcross==1) then
928 0 : ABI_FREE(ur_ae_bdgw)
929 0 : ABI_FREE(ur_ae_onsite_bdgw)
930 0 : ABI_FREE(ur_ps_onsite_bdgw)
931 : end if
932 : end if
933 : end do !spin
934 :
935 0 : ABI_FREE(igfftxg0)
936 0 : if (pawcross==1) then
937 0 : ABI_FREE(igfftfxg0)
938 : end if
939 : !
940 : ! Gather contributions from all the CPUs.
941 : !
942 : ! ===========================
943 : ! ==== Deallocate memory ====
944 : ! ===========================
945 0 : if (Psps%usepaw==1) then
946 : if (allocated(gwx_gfft)) then
947 : ABI_FREE(gwx_gfft)
948 : end if
949 0 : call pawcprj_free(Cprj_ksum)
950 0 : ABI_FREE(Cprj_ksum)
951 0 : if (allocated(Pwij_fft)) then
952 0 : call pawpwij_free(Pwij_fft)
953 0 : ABI_FREE(Pwij_fft)
954 : end if
955 : if (use_pawnhat==1) then
956 : ABI_FREE(nhat12)
957 : ABI_FREE(grnhat12)
958 : end if
959 0 : if (pawcross==1) then
960 0 : ABI_FREE(ur_ae_sum)
961 0 : ABI_FREE(ur_ae_onsite_sum)
962 0 : ABI_FREE(ur_ps_onsite_sum)
963 0 : ABI_FREE(ktabrf)
964 : end if
965 : end if
966 :
967 0 : ABI_FREE(wfr_sum)
968 0 : ABI_FREE(rhotwg_ki)
969 0 : ABI_FREE(vc_sqrt_qbz)
970 0 : ABI_FREE(ktabr)
971 : ! ABI_FREE(proc_distrb)
972 0 : if (plowan_compute<10) then
973 0 : ABI_FREE(coeffW_BZ)
974 : endif
975 :
976 :
977 0 : call timab(430,2,tsec) ! csigme (SigX)
978 :
979 : DBG_EXIT("COLL")
980 :
981 0 : end subroutine prep_calc_ucrpa
982 : !!***
983 :
984 : END MODULE m_prep_calc_ucrpa
985 : !!***
|
