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