Line data Source code
1 : !!****m* ABINIT/m_mlwfovlp_qp
2 : !! NAME
3 : !! m_mlwfovlp_qp
4 : !!
5 : !! FUNCTION
6 : !! Interpolate GW corrections with Wannier functions
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (DRH)
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 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_mlwfovlp_qp
23 :
24 : use defs_basis
25 : use defs_wannier90
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_hdr
30 : use m_dtset
31 : use m_dtfil
32 :
33 : use defs_abitypes, only : MPI_type
34 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
35 : use m_pawtab, only : pawtab_type
36 : use m_pawcprj, only : pawcprj_type, paw_overlap, pawcprj_getdim, pawcprj_alloc, pawcprj_free
37 : use m_pawrhoij, only : pawrhoij_type
38 : use m_numeric_tools, only : isordered
39 : use m_geometry, only : metric
40 : use m_crystal, only : crystal_t
41 : use m_kpts, only : listkk
42 : use m_bz_mesh, only : kmesh_t
43 : use m_ebands, only : ebands_t
44 : use m_qparticles, only : rdqps, rdgw
45 : use m_sort, only : sort_dp
46 :
47 : implicit none
48 :
49 : private
50 : !!***
51 :
52 : public :: mlwfovlp_qp
53 : !!***
54 :
55 : contains
56 : !!***
57 :
58 : !!****f* m_mlwfovlp_qp/mlwfovlp_qp
59 : !! NAME
60 : !! mlwfovlp_qp
61 : !!
62 : !! FUNCTION
63 : !! Routine which computes replaces DFT wave functions and eigenvalues with
64 : !! GW quasiparticle ones using previously computed qp wave functions in
65 : !! DFT bloch function representation for Wannier code (www.wannier.org f90 version).
66 : !!
67 : !! INPUTS
68 : !! dtset <type(dataset_type)>=all input variables for this dataset
69 : !! dtfil <type(datafiles_type)>=variables related to files
70 : !! mband=maximum number of bands
71 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
72 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
73 : !! mkmem =number of k points treated by this node.
74 : !! mpw=maximum dimensioned size of npw.
75 : !! natom=number of atoms in cell.
76 : !! nkpt=number of k points.
77 : !! npwarr(nkpt)=number of planewaves in basis at this k point
78 : !! nspden=number of spin-density components
79 : !! nsppol=1 for unpolarized, 2 for spin-polarized
80 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
81 : !! Hdr<Hdr_type>=The m_mlwfovlp_qp header.
82 : !! MPI_enreg=information about MPI parallelization
83 : !! Cprj_BZ(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
84 : !!
85 : !! OUTPUT
86 : !!
87 : !! SIDE EFFECTS
88 : !! cg(2,mcg)=planewave coefficients of wavefunctions
89 : !! replaced by quasiparticle wavefunctions
90 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues replaced by qp eigenvalues(hartree)
91 : !!
92 : !! NOTES
93 : !! Number of bands for wannier calculation must be identical to number used
94 : !! for gw calculation. Bands not wanted for wannier calculation must be
95 : !! excluded in exclude_band statement in wannier90.win file.
96 : !! Full plane-wave basis for DFT wavefunctions must be used in GW calculation,
97 : !! or inaccuracies may result.
98 : !! This is at best a beta version of this code, with little consistency
99 : !! checking, so the user must be very careful or the results may be invalid.
100 : !!
101 : !! SOURCE
102 :
103 1 : subroutine mlwfovlp_qp(cg,Cprj_BZ,dtset,dtfil,eigen,mband,mcg,mcprj,mkmem,mpw,natom,&
104 1 : & nkpt,npwarr,nspden,nsppol,ntypat,Hdr,Pawtab,rprimd,MPI_enreg)
105 :
106 : !Arguments ------------------------------------
107 : !scalars
108 : integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,nkpt,nspden,natom,ntypat
109 : integer,intent(in) :: nsppol
110 : type(dataset_type),intent(in) :: dtset
111 : type(datafiles_type),intent(in) :: dtfil
112 : type(Hdr_type),intent(in) :: Hdr
113 : type(MPI_type),intent(in) :: MPI_enreg
114 : type(pawcprj_type),target,intent(inout) :: Cprj_BZ(natom,mcprj)
115 : type(Pawtab_type),intent(in) :: Pawtab(ntypat*Dtset%usepaw)
116 : !arrays
117 : integer,intent(in) :: npwarr(nkpt)
118 : real(dp),intent(inout) :: cg(2,mcg)
119 : real(dp),intent(inout) :: eigen(mband*nkpt*nsppol)
120 : real(dp),intent(in) :: rprimd(3,3)
121 :
122 : !Local variables-------------------------------
123 : !scalars
124 : integer,parameter :: from_QPS_FILE=1,from_GW_FILE=2
125 : integer :: sppoldbl,timrev,bantot_ibz,ikibz,ikbz,dimrho
126 : integer :: iband,icg,icg_shift,ii,ipw,isppol,my_nspinor,nband_k,ord_iband
127 : integer :: nfftot,ikpt,irzkpt,npw_k,ikg
128 : integer :: nscf,nbsc,itimrev,band_index,nkibz,nkbz
129 : integer :: input !,jb_idx,ib_idx,ijpack, jband,
130 : integer :: nprocs,ios
131 : real(dp) :: TOL_SORT=tol12
132 : real(dp) :: dksqmax,ucvol !ortho_err,
133 : logical :: ltest,qpenek_is_ordered,g0w0_exists
134 : character(len=500) :: msg
135 : character(len=fnlen) :: gw_fname
136 1 : type(ebands_t) :: QP_bst
137 56 : type(crystal_t) :: Cryst
138 13 : type(kmesh_t) :: Kibz_mesh
139 1 : type(MPI_type) :: MPI_enreg_seq
140 : !arrays
141 2 : integer :: indkk(nkpt,6),my_ngfft(18)
142 2 : integer,allocatable :: npwarr_ibz(:),nband_ibz(:),ibz2bz(:,:),istwfk_ibz(:)
143 2 : integer,allocatable :: dimlmn(:),iord(:),nattyp_dum(:)
144 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) !,paw_ovlp(2)
145 1 : real(dp),allocatable :: qp_rhor(:,:),sorted_qpene(:)
146 2 : real(dp),allocatable :: kibz(:,:),wtk_ibz(:)
147 2 : real(dp),allocatable :: doccde_ibz(:),occfact_ibz(:),eigen_ibz(:)
148 1 : real(dp),allocatable :: igwene(:,:,:)
149 1 : complex(dp),allocatable :: m_ks_to_qp(:,:,:,:),m_ks_to_qp_BZ(:,:,:,:) !,ortho(:)
150 1 : complex(dp),allocatable :: m_tmp(:,:),cg_k(:,:),cg_qpk(:,:)
151 1 : type(Pawrhoij_type),allocatable :: prev_Pawrhoij(:)
152 : !type(pawcprj_type),pointer :: Cp1(:,:),Cp2(:,:)
153 : !************************************************************************
154 :
155 : ABI_UNUSED(mkmem)
156 :
157 : DBG_ENTER("COLL")
158 :
159 1 : write(msg,'(17a)')ch10,&
160 1 : ' mlwfovlp_qp: WARNING',ch10,&
161 1 : ' The input *_WFK file of DFT wavefunctions to be converted',ch10,&
162 1 : ' to GW quasiparticle wavefunctions MUST have been written in',ch10,&
163 1 : ' the run that produced the GW *_KSS file using kssform 3,',ch10,&
164 1 : ' the ONLY value of kssform permitted for GW Wannier functions.',ch10,&
165 1 : ' Otherwise, the *_QPS file needed here will be inconsistent,',ch10,&
166 1 : ' and the output quasiparticle wavefunctions will be garbage.',ch10,&
167 2 : ' No internal check that can verify this is presently possible.',ch10
168 1 : call wrtout(std_out,msg,'COLL')
169 :
170 : ! === Some features are not implemented yet ===
171 1 : ABI_CHECK(Dtset%nspinor==1,'nspinor==2 not implemented')
172 1 : ABI_CHECK(Dtset%nsppol==1,'nsppol==2 not implemented, check wannier90')
173 65 : ltest=ALL(Dtset%nband(1:Dtset%nkpt*Dtset%nsppol)==Dtset%nband(1))
174 1 : ABI_CHECK(ltest,'nband(:) should be constant')
175 : !
176 : ! MPI initialization
177 1 : nprocs=MPI_enreg%nproc_cell
178 :
179 1 : if (nprocs/=1) then
180 0 : ABI_ERROR("mlwfovlp_qp not programmed for parallel execution")
181 : end if
182 :
183 : ! Compute reciprocal space metric gmet for unit cell of disk wf
184 1 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
185 :
186 : ! Compute k points from gw irreducible set equivalent to full-zone wannier set
187 1 : sppoldbl=1 ; timrev=1 ; my_nspinor=max(1,Dtset%nspinor/MPI_enreg%nproc_spinor)
188 : call listkk(dksqmax,gmet,indkk,dtset%kptgw,dtset%kpt,dtset%nkptgw,nkpt,&
189 1 : & dtset%nsym,sppoldbl,dtset%symafm,dtset%symrel,timrev,xmpi_comm_self)
190 :
191 1 : if (dksqmax>tol8) then
192 : write(msg,'(5a)')&
193 0 : & 'Set of GW irreducible-zone kptgw in input file is inconsistent',ch10,&
194 0 : & 'with full-zone set being used for wannier90 setup.',ch10,&
195 0 : & 'Action: correct input data'
196 0 : ABI_ERROR(msg)
197 : end if
198 : !
199 : ! === Initialize object defining the Band strucuture ===
200 : ! * Initialize with KS results using IBZ indexing.
201 : ! * After rdqps, QP_bst will contain the QP amplitudes.
202 1 : nkibz=Dtset%nkptgw
203 3 : ABI_MALLOC(kibz,(3,nkibz))
204 3 : ABI_MALLOC(wtk_ibz,(nkibz))
205 34 : kibz=Dtset%kptgw(:,1:Dtset%nkptgw)
206 :
207 : ! MG: This part is needed to get the IBZ weight that will be reported
208 : ! on ab_out thus we should be consistent. Ideally Cryst should be
209 : ! one of the basic abinit objects and it should be passed to this routine.
210 :
211 : !different conventions are used in GW and abinit!!
212 1 : cryst = hdr%get_crystal(gw_timrev=timrev+1)
213 1 : call Kibz_mesh%init(Cryst,nkibz,kibz,Dtset%kptopt)
214 10 : wtk_ibz=Kibz_mesh%wt
215 1 : call cryst%free()
216 1 : call Kibz_mesh%free()
217 :
218 3 : ABI_MALLOC(ibz2bz,(nkibz,6))
219 : call listkk(dksqmax,gmet,ibz2bz,dtset%kpt,dtset%kptgw,nkpt,dtset%nkptgw,&
220 1 : & dtset%nsym,sppoldbl,dtset%symafm,dtset%symrel,timrev,xmpi_comm_self)
221 :
222 9 : ltest=ALL(ibz2bz(:,2)==1)
223 1 : ABI_CHECK(ltest,'Not able to found irreducible points in the BZ set!')
224 :
225 1 : if (dksqmax>tol8) then
226 : write(msg,'(5a)')&
227 0 : 'Set of GW irreducible-zone kptgw in input file is inconsistent',ch10,&
228 0 : 'with full-zone set being used for wannier90 setup.',ch10,&
229 0 : 'Action: correct input data'
230 0 : ABI_ERROR(msg)
231 : end if
232 :
233 3 : ABI_MALLOC(npwarr_ibz,(nkibz))
234 2 : ABI_MALLOC(istwfk_ibz,(nkibz))
235 3 : ABI_MALLOC(nband_ibz,(nkibz*nsppol))
236 :
237 2 : do isppol=1,nsppol
238 10 : do ikibz=1,nkibz
239 8 : ikbz=ibz2bz(ikibz+(sppoldbl-1)*(isppol-1)*nkibz,1)
240 8 : npwarr_ibz(ikibz)= npwarr(ikbz)
241 8 : istwfk_ibz(ikibz)=Dtset%istwfk(ikbz)
242 9 : nband_ibz(ikibz+(isppol-1)*nkibz)=Dtset%nband(ikbz+(isppol-1)*nkpt)
243 : end do
244 : end do
245 :
246 9 : bantot_ibz=SUM(nband_ibz)
247 3 : ABI_MALLOC(doccde_ibz,(bantot_ibz))
248 2 : ABI_MALLOC(eigen_ibz,(bantot_ibz))
249 2 : ABI_MALLOC(occfact_ibz,(bantot_ibz))
250 243 : doccde_ibz(:)=zero ; eigen_ibz(:)=zero ; occfact_ibz(:)=zero
251 :
252 : band_index=0
253 2 : do isppol=1,nsppol
254 10 : do ikibz=1,nkibz
255 8 : ikbz=ibz2bz(ikibz+(sppoldbl-1)*(isppol-1)*nkibz,1)
256 8 : nband_k=nband_ibz(ikibz+(isppol-1)*nkibz)
257 74 : ii=SUM(Dtset%nband(1:ikbz+(isppol-1)*nkpt))-nband_k
258 88 : eigen_ibz(band_index+1:band_index+nband_k)=eigen(ii+1:ii+nband_k)
259 9 : band_index=band_index+nband_k
260 : end do
261 : end do
262 :
263 : call QP_bst%init(bantot_ibz,Dtset%nelect,Dtset%ne_qFD,Dtset%nh_qFD,Dtset%ivalence,&
264 : doccde_ibz,eigen_ibz,istwfk_ibz,kibz,nband_ibz,&
265 : nkibz,npwarr_ibz,nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact_ibz,wtk_ibz,&
266 : dtset%cellcharge(1),dtset%kptopt,dtset%kptrlatt_orig,dtset%nshiftk_orig,dtset%shiftk_orig,&
267 1 : dtset%kptrlatt,dtset%nshiftk,dtset%shiftk)
268 :
269 1 : ABI_FREE(kibz)
270 1 : ABI_FREE(wtk_ibz)
271 1 : ABI_FREE(ibz2bz)
272 1 : ABI_FREE(npwarr_ibz)
273 1 : ABI_FREE(istwfk_ibz)
274 1 : ABI_FREE(nband_ibz)
275 1 : ABI_FREE(doccde_ibz)
276 1 : ABI_FREE(eigen_ibz)
277 1 : ABI_FREE(occfact_ibz)
278 :
279 : ! === Read in quasiparticle information ===
280 : ! * Initialize QP amplitudes with KS, QP_bst% presently contains KS energies.
281 : ! * If file not found return, everything has been already initialized with KS values
282 : ! Here qp_rhor is not needed thus dimrho=0
283 6 : ABI_MALLOC(m_ks_to_qp,(mband,mband,dtset%nkptgw,nsppol))
284 890 : m_ks_to_qp=czero
285 11 : do iband=1,mband
286 101 : m_ks_to_qp(iband,iband,:,:)=cone
287 : end do
288 :
289 : ! Fake MPI_type for rdqps
290 1 : call initmpi_seq(MPI_enreg_seq)
291 :
292 20 : my_ngfft=Dtset%ngfft; if (Dtset%usepaw==1.and.ALL(Dtset%ngfftdg(1:3)/=0)) my_ngfft=Dtset%ngfftdg
293 4 : nfftot=PRODUCT(my_ngfft(1:3)); dimrho=0
294 :
295 : ! Change gw_fname to read a GW file instead of the QPS file.
296 : ! TODO not so sure that wannier90 can handle G0W0 eigenvalues that are not ordered, though!
297 1 : gw_fname = "g0w0"
298 1 : g0w0_exists = .FALSE.
299 1 : inquire(file=gw_fname,iostat=ios,exist=g0w0_exists)
300 1 : if (ios/=0) then
301 0 : ABI_ERROR('File g0w0 exists but iostat returns nonzero value!')
302 : end if
303 :
304 1 : if (.not.g0w0_exists) then ! read QPS file (default behavior).
305 1 : input = from_QPS_FILE
306 3 : ABI_MALLOC(prev_Pawrhoij,(Cryst%natom*Dtset%usepaw))
307 2 : ABI_MALLOC(qp_rhor,(nfftot,nspden*dimrho))
308 :
309 : call rdqps(QP_bst,Dtfil%fnameabi_qps,Dtset%usepaw,Dtset%nspden,dimrho,nscf,&
310 1 : nfftot,my_ngfft,ucvol,Cryst,Pawtab,MPI_enreg_seq,nbsc,m_ks_to_qp,qp_rhor,prev_Pawrhoij)
311 :
312 1 : ABI_FREE(qp_rhor)
313 1 : ABI_FREE(prev_Pawrhoij)
314 :
315 : else
316 : ! Read GW file (m_ks_to_qp has been already set to 1, no extrapolation is performed)
317 0 : ABI_WARNING(' READING GW CORRECTIONS FROM FILE g0w0 !')
318 0 : input = from_GW_FILE
319 0 : ABI_MALLOC(igwene,(QP_bst%mband,QP_bst%nkpt,QP_bst%nsppol))
320 0 : call rdgw(QP_bst,gw_fname,igwene,extrapolate=.FALSE.)
321 0 : ABI_FREE(igwene)
322 : end if
323 :
324 : ! === Begin big loop over full-zone k points and spin (not implemented) ===
325 : ! * Wannier90 treats only a single spin, changes in wannier90 are needed
326 4 : ABI_MALLOC(cg_k,(mpw,mband))
327 3 : ABI_MALLOC(cg_qpk,(mpw,mband))
328 4 : ABI_MALLOC(m_tmp,(mband,mband))
329 :
330 2 : band_index=0 ; icg=0 ; ikg=0
331 2 : do isppol=1,nsppol
332 66 : do ikpt=1,nkpt
333 :
334 64 : irzkpt =indkk(ikpt+(sppoldbl-1)*(isppol-1)*nkpt,1)
335 64 : itimrev=indkk(ikpt+(sppoldbl-1)*(isppol-1)*nkpt,6)
336 64 : npw_k=npwarr(ikpt)
337 64 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
338 :
339 64 : if (nband_k/=mband) then
340 : write(msg,'(a,i0,7a)')&
341 0 : 'Number of bands for k point ',ikpt,' is inconsistent with number',ch10,&
342 0 : 'specified for wannier90 calculation',ch10,&
343 0 : 'Action: correct input so all band numbers are equal for GW',ch10,&
344 0 : 'and wannier90 datasets.'
345 0 : ABI_ERROR(msg)
346 : end if
347 :
348 : ! Load KS states for this kbz and spin
349 704 : do iband=1,nband_k
350 640 : icg_shift=npw_k*my_nspinor*(iband-1)+icg
351 121214 : do ipw=1,npw_k
352 121150 : cg_k(ipw,iband)=DCMPLX(cg(1,ipw+icg_shift),cg(2,ipw+icg_shift))
353 : end do
354 : end do
355 :
356 : ! If time reversal is used for relating ikpt to irzkpt, then multiply by
357 : ! the complex conjugate of the lda-to-qp transformation matrix
358 64 : if (itimrev==0) then
359 5328 : m_tmp(:,:)=m_ks_to_qp(:,:,irzkpt,isppol)
360 16 : else if (itimrev==1) then
361 1776 : m_tmp(:,:)=conjg(m_ks_to_qp(:,:,irzkpt,isppol))
362 : else
363 0 : write(msg,'(2(a,i0))')'Invalid indkk(ikpt,6) ',itimrev,'from routine listkk for k-point ',ikpt
364 0 : ABI_BUG(msg)
365 : end if
366 :
367 64 : call ZGEMM('N','N',npw_k,mband,mband,cone,cg_k,mpw,m_tmp,mband,czero,cg_qpk,mpw)
368 :
369 : ! === Orthonormality test ===
370 : ! * nband >= maxval(bndgw) for this to pass, but may be less than nband used in GW.
371 : ! * Unfortunately, does not test WFK and QPS consistency.
372 : !allocate(ortho(nband_k*(nband_k+1)/2))
373 : !ortho=czero; ijpack=0
374 : !do jband=1,nband_k
375 : ! jb_idx=band_index+jband
376 : ! if (dtset%usepaw==1) Cp2 => Cprj_BZ(:,jband:jband+(my_nspinor-1))
377 : ! do iband=1,jband
378 : ! ib_idx=band_index+iband
379 : ! ijpack=ijpack+1
380 : ! ortho(ijpack)=sum(conjg(cg_qpk(1:npw_k,iband))*cg_qpk(1:npw_k,jband))
381 : ! if (dtset%usepaw==1) then
382 : ! Cp1 => Cprj_BZ(:,iband:iband+(my_nspinor-1))
383 : ! paw_ovlp = paw_overlap(Cp2,Cp1,Cryst%typat,Pawtab)
384 : ! ortho(ijpack) = ortho(ijpack) + CMPLX(paw_ovlp(1),paw_ovlp(2))
385 : ! end if
386 : ! if (jband==iband) ortho(ijpack)=ortho(ijpack)-cone
387 : ! end do
388 : !end do
389 : !ortho_err=maxval(abs(ortho))
390 :
391 : !write(std_out,*)' drh - mlwfovlp_qp: ikpt,ortho_err',ikpt,ortho_err
392 : !if (ortho_err>tol6) then
393 : ! write(msg, '(3a,i4,a,i6,a,1p,e8.1,3a)' )&
394 : !& ' orthonormality error for quasiparticle wave functions.',ch10,&
395 : !& ' spin=',isppol,' k point=',ikpt,' ortho_err=',ortho_err,' >1E-6',ch10,&
396 : !& ' Action: Be sure input nband>=maxval(bndgw)'
397 : ! ABI_ERROR(msg)
398 : !end if
399 : !deallocate(ortho)
400 :
401 : ! Replace lda wave functions and eigenvalues with quasiparticle ones.
402 64 : qpenek_is_ordered = isordered(nband_k,QP_bst%eig(:,irzkpt,isppol),">",TOL_SORT)
403 :
404 64 : if (input==from_QPS_FILE .and. .not.qpenek_is_ordered) then
405 : write(msg,'(3a)')&
406 29 : " QP energies read from QPS file are not ordered, likely nband_k>nbdgw. ",ch10,&
407 58 : " Change nband in the input file so that it equals the number of GW states calculated"
408 29 : ABI_WARNING(msg)
409 : end if
410 :
411 : if ( .TRUE. ) then
412 704 : do iband=1,nband_k
413 640 : icg_shift=npw_k*my_nspinor*(iband-1)+icg
414 640 : eigen(iband+band_index)=QP_bst%eig(iband,irzkpt,isppol)
415 121214 : do ipw=1,npw_k
416 120510 : cg(1,ipw+icg_shift)= real(cg_qpk(ipw,iband))
417 121150 : cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,iband))
418 : end do
419 : end do
420 : else
421 : ! FIXME There's a problem in twannier90 since nband_k > nbdgw and therefore we also read KS states from the QPS file!
422 : ! Automatic test has to be changed!
423 : write(msg,'(2a,3f8.4,3a)')ch10,&
424 : "QP energies at k-point ",QP_bst%kptns(:,irzkpt)," are not sorted in ascending numerical order!",ch10,&
425 : "Performing reordering of energies and wavefunctions to be written on the final WKF file."
426 : ABI_ERROR(msg)
427 : !write(std_out,*)"eig",(QP_bst%eig(ii,irzkpt,isppol),ii=1,nband_k)
428 : ABI_MALLOC(sorted_qpene,(nband_k))
429 : ABI_MALLOC(iord,(nband_k))
430 : sorted_qpene = QP_bst%eig(1:nband_k,irzkpt,isppol)
431 : iord = (/(ii, ii=1,nband_k)/)
432 :
433 : call sort_dp(nband_k,sorted_qpene,iord,TOL_SORT)
434 : do ii=1,nband_k
435 : write(std_out,*)"%eig, sorted_qpene, iord",QP_bst%eig(ii,irzkpt,isppol)*Ha_eV,sorted_qpene(ii)*Ha_eV,iord(ii)
436 : end do
437 :
438 : do iband=1,nband_k
439 : ord_iband = iord(iband)
440 : icg_shift=npw_k*my_nspinor*(iband-1)+icg
441 : !eigen(iband+band_index)=QP_bst%eig(iband,irzkpt,isppol)
442 : eigen(iband+band_index)=QP_bst%eig(ord_iband,irzkpt,isppol)
443 : do ipw=1,npw_k
444 : !cg(1,ipw+icg_shift)= real(cg_qpk(ipw,iband))
445 : !cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,iband))
446 : cg(1,ipw+icg_shift)= real(cg_qpk(ipw,ord_iband))
447 : cg(2,ipw+icg_shift)=aimag(cg_qpk(ipw,ord_iband))
448 : end do
449 : end do
450 : ABI_FREE(sorted_qpene)
451 : ABI_FREE(iord)
452 : end if
453 :
454 64 : band_index=band_index+nband_k
455 64 : icg=icg+npw_k*my_nspinor*nband_k
456 1 : ikg=ikg+npw_k
457 : end do !ikpt
458 : end do !isppol
459 :
460 1 : ABI_FREE(cg_k)
461 1 : ABI_FREE(cg_qpk)
462 1 : ABI_FREE(m_tmp)
463 :
464 : ! === If PAW, update projections in BZ ===
465 : ! * Since I am lazy and here I do not care about memory, I just reconstruct m_ks_to_qp in the BZ.
466 : ! * update_cprj will take care of updating the PAW projections to get <p_lmn|QP_{nks]>
467 : ! This allows some CPU saving, no need to call ctocprj.
468 : ! FIXME this part should be tested, automatic test to be provided
469 1 : if (Dtset%usepaw==1) then
470 0 : ABI_MALLOC(dimlmn,(natom))
471 0 : call pawcprj_getdim(dimlmn,dtset%natom,nattyp_dum,ntypat,Dtset%typat,pawtab,'R')
472 :
473 0 : nkbz=nkpt
474 0 : ABI_MALLOC(m_ks_to_qp_BZ,(mband,mband,nkbz,nsppol))
475 0 : do isppol=1,nsppol
476 0 : do ikbz=1,nkbz
477 0 : ikibz =indkk(ikibz+(sppoldbl-1)*(isppol-1)*nkbz,1)
478 0 : itimrev=indkk(ikibz+(sppoldbl-1)*(isppol-1)*nkbz,6)
479 0 : select case (itimrev)
480 : case (0)
481 0 : m_ks_to_qp_BZ(:,:,ikbz,isppol)=m_ks_to_qp(:,:,ikibz,isppol)
482 : case (1)
483 0 : m_ks_to_qp_BZ(:,:,ikbz,isppol)=CONJG(m_ks_to_qp(:,:,ikibz,isppol))
484 : case default
485 0 : write(msg,'(a,i3)')"Wrong itimrev= ",itimrev
486 0 : ABI_BUG(msg)
487 : end select
488 : end do
489 : end do
490 :
491 0 : call update_cprj(natom,nkbz,mband,nsppol,my_nspinor,m_ks_to_qp_BZ,dimlmn,Cprj_BZ)
492 0 : ABI_FREE(dimlmn)
493 0 : ABI_FREE(m_ks_to_qp_BZ)
494 : end if !PAW
495 :
496 1 : write(msg,'(6a)')ch10,&
497 1 : ' mlwfovlp_qp: Input KS wavefuctions have been converted',ch10,&
498 1 : ' to GW quasiparticle wavefunctions for maximally localized wannier',ch10,&
499 2 : ' function construction by wannier90.'
500 1 : call wrtout(ab_out,msg,'COLL')
501 1 : call wrtout(std_out,msg,'COLL')
502 :
503 1 : ABI_FREE(m_ks_to_qp)
504 1 : call QP_bst%free()
505 1 : call destroy_mpi_enreg(MPI_enreg_seq)
506 :
507 : DBG_EXIT("COLL")
508 :
509 1 : end subroutine mlwfovlp_qp
510 : !!***
511 :
512 : !!****f* ABINIT/update_cprj
513 : !! NAME
514 : !! update_cprj
515 : !!
516 : !! FUNCTION
517 : !! Update the matrix elements of the PAW projectors in case of self-consistent GW.
518 : !!
519 : !! INPUTS
520 : !! dimlmn(natom)=number of (l,m,n) components for each atom (only for PAW)
521 : !! nkibz=number of k-points
522 : !! nsppol=number of spin
523 : !! nbnds=number of bands in the present GW calculation
524 : !! m_ks_to_qp(nbnds,nbnds,nkibz,nsppol)= expansion of the QP amplitudes in terms of KS wavefunctions
525 : !! natom=number of atomd in unit cell
526 : !!
527 : !! OUTPUT
528 : !! Cprj_ibz(natom,nspinor*nkibz*nbnds*nsppol) <type(pawcprj_type)>=projected wave functions
529 : !! <Proj_i|Cnk> with all NL projectors. On exit, it contains the projections onto the
530 : !! QP amplitudes.
531 : !!
532 : !! TODO
533 : !! To be moved to cprj_utils, although here we use complex variables.
534 : !!
535 : !! SOURCE
536 : !!
537 :
538 0 : subroutine update_cprj(natom,nkibz,nbnds,nsppol,nspinor,m_ks_to_qp,dimlmn,Cprj_ibz)
539 :
540 : !Arguments ------------------------------------
541 : !scalars
542 : integer,intent(in) :: natom,nbnds,nkibz,nsppol,nspinor
543 : !arrays
544 : integer,intent(in) :: dimlmn(natom)
545 : complex(dp),intent(in) :: m_ks_to_qp(nbnds,nbnds,nkibz,nsppol)
546 : type(pawcprj_type),intent(inout) :: Cprj_ibz(natom,nspinor*nbnds*nkibz*nsppol)
547 :
548 : !Local variables-------------------------------
549 : !scalars
550 : integer :: iat,ib,ik,is,shift,indx_kibz,ilmn,nlmn,ispinor,ibsp,spad,ibdx
551 : !arrays
552 0 : real(dp),allocatable :: re_p(:),im_p(:),vect(:,:),umat(:,:,:)
553 : type(pawcprj_type),allocatable :: Cprj_ks(:,:)
554 : !************************************************************************
555 :
556 : DBG_ENTER("COLL")
557 :
558 0 : ABI_MALLOC(Cprj_ks,(natom,nspinor*nbnds))
559 0 : call pawcprj_alloc(Cprj_ks,0,dimlmn)
560 :
561 0 : ABI_MALLOC(re_p,(nbnds))
562 0 : ABI_MALLOC(im_p,(nbnds))
563 0 : ABI_MALLOC(vect,(2,nbnds))
564 0 : ABI_MALLOC(umat,(2,nbnds,nbnds))
565 : !
566 : ! $ \Psi^{QP}_{r,b} = \sum_n \Psi^{KS}_{r,n} M_{n,b} $
567 : !
568 : ! therefore the updated PAW projections are given by:
569 : !
570 : ! $ \<\tprj_j|\Psi^{QP}_a\> = sum_b M_{b,a} <\tprj_j|\Psi^{KS}_b\> $.
571 : !
572 0 : do is=1,nsppol
573 0 : do ik=1,nkibz
574 :
575 0 : shift=nspinor*nbnds*nkibz*(is-1)
576 0 : indx_kibz=nspinor*nbnds*(ik-1)+shift
577 0 : ibsp=0
578 0 : do ib=1,nbnds
579 0 : do ispinor=1,nspinor
580 0 : ibsp=ibsp+1
581 0 : do iat=1,natom
582 0 : Cprj_ks(iat,ibsp)%cp(:,:)=Cprj_ibz(iat,indx_kibz+ibsp)%cp(:,:)
583 : end do
584 : end do
585 : end do
586 :
587 0 : umat(1,:,:)=TRANSPOSE( REAL (m_ks_to_qp(:,:,ik,is)) )
588 0 : umat(2,:,:)=TRANSPOSE( AIMAG(m_ks_to_qp(:,:,ik,is)) )
589 :
590 0 : do iat=1,natom
591 0 : nlmn=dimlmn(iat)
592 0 : do ilmn=1,nlmn
593 :
594 0 : do ispinor=1,nspinor
595 : ! * Retrieve projections for this spinor component, at fixed atom and ilmn.
596 0 : spad=(ispinor-1)
597 0 : ibdx=0
598 0 : do ib=1,nbnds*nspinor,nspinor
599 0 : ibdx=ibdx+1
600 0 : vect(1,ibdx)=Cprj_ks(iat,ib+spad)%cp(1,ilmn)
601 0 : vect(2,ibdx)=Cprj_ks(iat,ib+spad)%cp(2,ilmn)
602 : end do
603 :
604 : re_p(:)= &
605 : & MATMUL(umat(1,:,:),vect(1,:)) &
606 0 : & -MATMUL(umat(2,:,:),vect(2,:))
607 :
608 : im_p(:)= &
609 : & MATMUL(umat(1,:,:),vect(2,:)) &
610 0 : & +MATMUL(umat(2,:,:),vect(1,:))
611 :
612 : ! === Save values ===
613 0 : ibdx=0
614 0 : do ib=1,nbnds*nspinor,nspinor
615 0 : ibdx=ibdx+1
616 0 : Cprj_ibz(iat,indx_kibz+spad+ib)%cp(1,ilmn)=re_p(ibdx)
617 0 : Cprj_ibz(iat,indx_kibz+spad+ib)%cp(2,ilmn)=im_p(ibdx)
618 : end do
619 : end do !ispinor
620 :
621 : end do !ilmn
622 : end do !iat
623 :
624 : end do !ik
625 : end do !is
626 :
627 0 : ABI_FREE(re_p)
628 0 : ABI_FREE(im_p)
629 0 : ABI_FREE(vect)
630 0 : ABI_FREE(umat)
631 :
632 0 : call pawcprj_free(Cprj_ks)
633 0 : ABI_FREE(Cprj_ks)
634 :
635 : DBG_EXIT("COLL")
636 :
637 0 : end subroutine update_cprj
638 : !!***
639 :
640 0 : end module m_mlwfovlp_qp
641 : !!***
|
