Branch data 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:Gygi1986]] 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 : : !!***
|
