Line data Source code
1 : !!****m* ABINIT/m_qparticles
2 : !! NAME
3 : !! m_qparticles
4 : !!
5 : !! FUNCTION
6 : !! This module contains tools for the IO of the QP file and other procedures
7 : !! related to the calculation of the quasiparticle amplitudes represented in terms
8 : !! of KS states.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2008-2026 ABINIT group (FB, MG)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : MODULE m_qparticles
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_hdr
29 : use m_errors
30 : use m_nctk
31 :
32 : use defs_datatypes, only : pseudopotential_type
33 : use defs_abitypes, only : MPI_type
34 : use m_io_tools, only : open_file, file_exists, isncfile
35 : use m_fstrings, only : int2char10, itoa, sjoin
36 : use m_numeric_tools, only : linfit, c2r, set2unit, interpol3d_0d, rhophi
37 : use m_gwdefs, only : sigparams_t
38 : use m_crystal, only : crystal_t
39 : use m_bz_mesh, only : kmesh_t
40 : use m_ebands, only : ebands_t
41 : use m_sigma, only : sigma_t
42 : use m_pawtab, only : pawtab_type
43 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_io, pawrhoij_inquire_dim
44 : use m_fourier_interpol,only : fourier_interpol
45 :
46 : implicit none
47 :
48 : private
49 :
50 : public :: wrqps ! Write a QPS file.
51 : public :: rdqps ! Read a QPS file.
52 : public :: show_QP ! Report the components of a QP amplitude in terms of KS eigenstates.
53 : public :: rdgw ! Read GW corrections from an external file.
54 : public :: updt_m_ks_to_qp ! Update the matrix of unitary transformation from lda to qp states.
55 :
56 : CONTAINS !=======================================================================================
57 : !!***
58 :
59 : !!****f* m_qparticles/wrqps
60 : !! NAME
61 : !! wrqps
62 : !!
63 : !! FUNCTION
64 : !! Write the _QPS file containing information on the quasi-particles energies and wavefunctions.
65 : !!
66 : !! INPUTS
67 : !! fname=The name of the file
68 : !! Sigp<sigparams_t>=Parameters characterizing the self-energy calculation.
69 : !! %nsppol=1 for unpolarized, 2 for spin-polarized
70 : !! %nbnds=number of bands used for sigma
71 : !! Sr<sigma_t>=Structure containing the results of the sigma run.
72 : !! %en_qp_diago(nbnds,nibz,nsppol)= NEW quasi-particle energies
73 : !! %eigvec_qp(nbnds,nbnds,nibz,nsppol)= NEW QP amplitudes in the KS basis set
74 : !! obtained by diagonalizing H0 + Herm(Sigma).
75 : !! m_ks_to_qp(nbnds,nbnds,nibz,nsppol)= expansion of the OLD QP amplitudes in terms of KS wavefunctions
76 : !! Kmesh<kmesh_t>=information on the k-point sampling.
77 : !! %nibz=number of irreducible k-points
78 : !! %ibz(3,kibz)=reduced coordinates of the irreducible k-points
79 : !! nfftot=Total number of FFT points for density
80 : !! ngfftf(18)=Info on the FFT mesh for the density.
81 : !! nscf=Number of self consistent cycles performed
82 : !! nspden=number of spin-density components
83 : !! Cryst<crystal_t>=Structure defining the crystal structure.
84 : !! Psps<type(pseudopotential_type)>=variables related to pseudopotentials.
85 : !! Pawrhoij(Cryst%natom*Psps%usepaw)<type(pawrhoij_type)>= rhoij datastructure.
86 : !! BSt<ebands_t>=Structure containing the band structure energies (only used is nscf==-1)
87 : !!
88 : !! OUTPUT
89 : !! Only writing
90 : !!
91 : !! NOTES
92 : !! Old QPS fileformat:
93 : !! |
94 : !! | No. of QPSCF cycles already performed.
95 : !! | No. of k-points in the IBZ.
96 : !! | Total number of bands used to construct the Green's function (nbnds)
97 : !! | nsppol
98 : !! | For each spin and k-point in the IBZ:
99 : !! | Reduced coordinates of the k-point.
100 : !! | for each band:
101 : !! | QP energies obtained by diagonalizing the QPSCGW Hamiltonian.
102 : !! | <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>$, ib=1,nbnds
103 : !! | FFT dimensions of the fine grid
104 : !! | QP density in real space.
105 : !!
106 : !! SOURCE
107 :
108 55 : subroutine wrqps(fname,Sigp,Cryst,Kmesh,Psps,Pawtab,Pawrhoij,nspden,nscf,nfftot,ngfftf,Sr,Bst,m_ks_to_qp,rho_qp)
109 :
110 : !Arguments ------------------------------------
111 : !scalars
112 : integer,intent(in) :: nfftot,nscf,nspden
113 : character(len=*),intent(in) :: fname
114 : type(kmesh_t),intent(in) :: Kmesh
115 : type(ebands_t),intent(in) :: BSt
116 : type(sigparams_t),intent(in) :: Sigp
117 : type(sigma_t),intent(in) :: Sr
118 : type(crystal_t),intent(in) :: Cryst
119 : type(Pseudopotential_type),intent(in) :: Psps
120 : !arrays
121 : integer,intent(in) :: ngfftf(18)
122 : real(dp),intent(in) :: rho_qp(nfftot,nspden)
123 : complex(dp),intent(in) :: m_ks_to_qp(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
124 : type(Pawrhoij_type),intent(inout) :: Pawrhoij(Cryst%natom*Psps%usepaw)
125 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
126 :
127 : !Local variables-------------------------------
128 : !scalars
129 : integer :: ib,ik,is,unqps,iatom,itypat
130 : character(len=500) :: msg
131 : !arrays
132 55 : integer,allocatable :: nlmn_type(:)
133 55 : complex(dp),allocatable :: mtmp(:,:)
134 : ! *************************************************************************
135 :
136 : DBG_ENTER("COLL")
137 :
138 55 : if (nscf >= 0) then
139 55 : write(msg,'(3a)')ch10,' writing QP data on file : ',TRIM(fname)
140 165 : call wrtout([std_out, ab_out], msg)
141 : end if
142 :
143 55 : if (open_file(fname,msg,newunit=unqps,form='formatted',status='unknown') /= 0) then
144 0 : ABI_ERROR(msg)
145 : end if
146 :
147 55 : write(unqps,*)nscf+1
148 55 : write(unqps,*)Kmesh%nibz
149 55 : write(unqps,*)Sigp%nbnds
150 55 : write(unqps,*)Sigp%nsppol
151 :
152 220 : ABI_MALLOC(mtmp,(Sigp%nbnds,Sigp%nbnds))
153 :
154 55 : if (nscf>=0) then
155 : ! Write the new m_ks_to_qp on file.
156 112 : do is=1,Sigp%nsppol
157 463 : do ik=1,Kmesh%nibz
158 351 : write(unqps,*)Kmesh%ibz(:,ik)
159 4688 : do ib=1,Sigp%nbnds
160 4280 : write(unqps,*)Sr%en_qp_diago(ib,ik,is)
161 4631 : write(unqps,*)m_ks_to_qp(:,ib,ik,is)
162 : end do
163 : end do
164 : end do
165 0 : else if (nscf==-1) then
166 : ! Write fake QPS file with KS band structure (Mainly used for G0W)
167 0 : call set2unit(mtmp)
168 0 : do is=1,Sigp%nsppol
169 0 : do ik=1,Kmesh%nibz
170 0 : write(unqps,*)Kmesh%ibz(:,ik)
171 0 : do ib=1,Sigp%nbnds
172 0 : write(unqps,*)BSt%eig(ib,ik,is)
173 0 : write(unqps,*)mtmp(:,ib)
174 : end do
175 : end do
176 : end do
177 : else
178 0 : ABI_ERROR(sjoin("Wrong nscf ",itoa(nscf)))
179 : end if
180 :
181 55 : ABI_FREE(mtmp)
182 :
183 1017104 : write(msg,'(a,f9.4)')' (wrqps) planewave contribution to nelect: ',SUM(rho_qp(:,1))*Cryst%ucvol/nfftot
184 55 : call wrtout(std_out,msg)
185 55 : if (nspden == 4) then
186 0 : write(msg,'(a,3f9.4)')' mx, my, mz: ',&
187 0 : SUM(rho_qp(:,2))*Cryst%ucvol/nfftot,SUM(rho_qp(:,3))*Cryst%ucvol/nfftot,SUM(rho_qp(:,4))*Cryst%ucvol/nfftot
188 0 : call wrtout(std_out,msg)
189 : end if
190 :
191 : ! Write FFT dimensions and QP density
192 55 : write(unqps,*)ngfftf(1:3)
193 55 : write(unqps,*)rho_qp(:,:)
194 :
195 55 : if (Psps%usepaw==1) then
196 : ! Write QP rhoij to be used for on-site density mixing.
197 0 : ABI_MALLOC(nlmn_type,(Cryst%ntypat))
198 0 : do itypat=1,Cryst%ntypat
199 0 : nlmn_type(itypat)=Pawtab(itypat)%lmn_size
200 : end do
201 :
202 0 : write(unqps,*) Cryst%natom, Cryst%ntypat
203 0 : write(unqps,*) (Cryst%typat(iatom), iatom=1,Cryst%natom)
204 0 : write(unqps,*) (nlmn_type(itypat), itypat=1,Cryst%ntypat)
205 0 : write(unqps,*) Pawrhoij(1)%nsppol, Pawrhoij(1)%nspden
206 :
207 : call pawrhoij_io(pawrhoij,unqps,Sigp%nsppol,Sigp%nspinor,nspden,nlmn_type,Cryst%typat,&
208 0 : & HDR_LATEST_HEADFORM,"Write",form="formatted")
209 0 : ABI_FREE(nlmn_type)
210 : end if
211 :
212 55 : close(unqps)
213 :
214 : DBG_EXIT("COLL")
215 :
216 55 : end subroutine wrqps
217 : !!***
218 :
219 : !----------------------------------------------------------------------
220 :
221 : !!****f* m_qparticles/rdqps
222 : !! NAME
223 : !! rdqps
224 : !!
225 : !! FUNCTION
226 : !! Read a _QPS file containing the QP energies of the previous iteration, the coefficients
227 : !! defining the QP amplitudes in terms of the KS basis set and the QP density for mixing.
228 : !!
229 : !! INPUTS
230 : !! nfftot=Total number of FFT points for density
231 : !! ngfftf(18)=Info on the FFT mesh for the density.
232 : !! nspden=Number of SPin-DENsity components.
233 : !! usepaw=1 if we are using PAW.
234 : !! fname=Name of the file
235 : !! dimrho=1 if density has to be read, 0 otherwise
236 : !! BSt<ebands_t>=Structure containing the initial band structure.
237 : !! ucvol=Volume of the unit cell
238 : !!
239 : !! OUTPUT
240 : !! nbsc=number of bands used to describe the QP amplitudes
241 : !! nscf=number of iterations that have been performed (==0 if we start from a KS calculation)
242 : !! m_ks_to_qp(mband,mband,nibz,nsppol)=matrix giving the decomposition of the QP
243 : !! wavefunction in the mainfold generated by the KS wavefunctions
244 : !! (i.e. $ m_ks_to_qp(ib,jb,k,s) := <\psi_{ib,k,s}^{KS}|\psi_{jb,k,s}^{QP}>$
245 : !! rhor_out(nfftot,nspden)=quasiparticle density
246 : !!
247 : !! SIDE EFFECTS
248 : !! BSt<ebands_t>=Structure containing the initial band structure.
249 : !! %en_qp(mband,nkpt,nsppol)=QP energies at iteration nscf
250 : !!
251 : !! TODO
252 : !! The value of nspden is not reported in the QPS file thus we have a possible undetected error.
253 : !!
254 : !! SOURCE
255 :
256 126 : subroutine rdqps(BSt,fname,usepaw,nspden,dimrho,nscf,&
257 126 : & nfftot,ngfftf,ucvol,Cryst,Pawtab,MPI_enreg,nbsc,m_ks_to_qp,rhor_out,Pawrhoij)
258 :
259 : !Arguments ------------------------------------
260 : !scalars
261 : integer,intent(in) :: nfftot,nspden,usepaw,dimrho
262 : integer,intent(out) :: nbsc,nscf
263 : real(dp),intent(in) :: ucvol
264 : character(len=*),intent(in) :: fname
265 : type(crystal_t),intent(in) :: Cryst
266 : type(ebands_t),intent(inout) :: BSt
267 : type(MPI_type),intent(inout) :: MPI_enreg
268 : !arrays
269 : integer,intent(in) :: ngfftf(18)
270 : real(dp),intent(out) :: rhor_out(nfftot,nspden*dimrho)
271 : complex(dp),intent(out) :: m_ks_to_qp(BSt%mband,BSt%mband,BSt%nkpt,BSt%nsppol)
272 : type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*usepaw)
273 : type(Pawrhoij_type),intent(inout) :: Pawrhoij(Cryst%natom*usepaw)
274 :
275 : !Local variables-------------------------------
276 : !scalars
277 : integer,parameter :: master=0
278 : integer :: ib,ii,ik,isppol,nbandR,nkibzR,nsppolR,unqps,my_rank,ispden
279 : integer :: ifft,n1,n2,n3,ir1,ir2,ir3,ios
280 : integer :: cplex_fft,optin,optout,nfft_found
281 : integer :: iatom,natomR,nspdenR,ntypatR,itypat
282 : real(dp) :: uerr,nelect_qps,ratio
283 : logical,parameter :: use_FFT_interpolation=.TRUE.
284 : logical :: ltest
285 : character(len=500) :: msg
286 : !arrays
287 : integer :: ngfft_found(18)
288 126 : integer,allocatable :: nlmn_type(:),typatR(:)
289 : real(dp) :: kibz(3),rr(3),rhogdum(1,1)
290 126 : real(dp),allocatable :: en_tmp(:), rhor_tmp(:,:)
291 126 : complex(dp),allocatable :: mtmp(:,:),utest(:,:)
292 : ! *************************************************************************
293 :
294 : DBG_ENTER("COLL")
295 :
296 837 : ABI_CHECK(ALL(BSt%nband==BSt%nband(1)), "No. of bands must be constant")
297 126 : ABI_CHECK(dimrho==0.or.dimrho==1, 'dimrho must be 0 or 1')
298 :
299 : ! This does not work in parallel !!?
300 : !% my_rank = xmpi_comm_rank(MPI_enreg%spaceComm)
301 126 : my_rank = MPI_enreg%me_kpt
302 :
303 : ! Check whether file exists or not.
304 126 : write(msg,'(5a)')ch10,&
305 126 : ' rdqps: reading QP wavefunctions of the previous step ',ch10,&
306 252 : ' looking for file ',TRIM(fname)
307 378 : call wrtout([std_out, ab_out], msg)
308 :
309 126 : if (.not.file_exists(fname)) then
310 64 : write(msg,'(2a)')' file not found, 1st iteration initialized with KS eigenelements ',ch10
311 192 : call wrtout([std_out, ab_out], msg)
312 128 : nscf=0; RETURN
313 : end if
314 :
315 62 : if (.not.isncfile(fname)) then
316 62 : if (open_file(fname,msg,newunit=unqps,form='formatted',status='unknown') /= 0) then
317 0 : ABI_ERROR(msg)
318 : end if
319 :
320 : ! TODO the _QPS file should contain additional information
321 62 : read(unqps,*)nscf
322 62 : write(msg,'(a,i4,a)')' Number of iteration(s) already performed: ',nscf,ch10
323 186 : call wrtout([std_out, ab_out], msg)
324 :
325 62 : read(unqps,*)nkibzR
326 62 : if (nkibzR/=BSt%nkpt) then
327 0 : write(msg,'(2(a,i0))')'Wrong number of k-points; Expected: ',BSt%nkpt,', Found: ',nkibzR
328 0 : ABI_ERROR(msg)
329 : end if
330 :
331 62 : read(unqps,*)nbandR
332 62 : nbsc=MIN(nbandR,BSt%mband)
333 :
334 62 : if (nbsc/=BSt%mband) then
335 : write(msg,'(3a,i4,a,i4)')&
336 0 : 'QPS file contains less bands than that used in the present calculation ',ch10,&
337 0 : 'Required: ',BSt%mband,', Found: ',nbandR
338 0 : ABI_WARNING(msg)
339 : end if
340 :
341 62 : if (nbsc/=nbandR) then
342 : write(msg,'(3a,i0,a)')&
343 0 : 'The QPS file contains more bands than that used in the present calculation ',ch10,&
344 0 : 'only the first ',nbandR,' bands will be read'
345 0 : ABI_COMMENT(msg)
346 : end if
347 :
348 248 : ABI_MALLOC(mtmp,(nbandR,nbandR))
349 186 : ABI_MALLOC(en_tmp,(nbandR))
350 62 : read(unqps,*)nsppolR
351 :
352 62 : ABI_CHECK_IEQ(nsppolR, BSt%nsppol, "QPS generated with different nsppol")
353 :
354 : ! Read energies and transformation for each k-point and spin.
355 : ! TODO: The format of the QPS file must be standardized !
356 : ! For example we might add the occupation numbers.
357 124 : do isppol=1,BSt%nsppol
358 456 : do ik=1,BSt%nkpt
359 332 : read(unqps,*)kibz(:)
360 332 : write(msg,'(a,i5,a,3(f6.3,1x),4x,a,i2)')' Reading ik ',ik,') k = ',kibz(:),' is = ',isppol
361 332 : call wrtout(std_out,msg)
362 1328 : ltest=(ALL(ABS(kibz(:)-BSt%kptns(:,ik))<0.001))
363 332 : ABI_CHECK(ltest,'Wrong k-point read')
364 4396 : do ib=1,nbandR
365 4064 : read(unqps,*)en_tmp(ib)
366 4396 : read(unqps,*)mtmp(:,ib)
367 : end do
368 :
369 : ! Store transformation and update energies.
370 55980 : m_ks_to_qp(1:nbsc,1:nbsc,ik,isppol)=mtmp(1:nbsc,1:nbsc)
371 4396 : BSt%eig(1:nbsc,ik,isppol)=en_tmp(1:nbsc)
372 :
373 : ! Check if matrix is unitary.
374 1328 : ABI_MALLOC(utest,(nbsc,nbsc))
375 55980 : utest(:,:) = TRANSPOSE(mtmp(1:nbsc,1:nbsc)) !this is just for the buggy gfortran
376 846700 : utest(:,:) = MATMUL(CONJG(utest),mtmp(1:nbsc,1:nbsc))
377 4396 : do ii=1,nbsc
378 4396 : utest(ii,ii)=utest(ii,ii)-one
379 : end do
380 55980 : uerr=MAXVAL(ABS(utest))
381 332 : if (uerr>tol6) then
382 0 : write(msg,'(a,es16.8)')' KS -> QP matrix is not unitary, MAX error = ',uerr
383 0 : ABI_WARNING(msg)
384 : end if
385 394 : ABI_FREE(utest)
386 : end do !ik
387 : end do !isppol
388 :
389 62 : ABI_FREE(mtmp)
390 62 : ABI_FREE(en_tmp)
391 :
392 : ! Read the QP density.
393 : ! The two FFT grids might differ. In case perform an FFT interpolation to have rhor on the input mesh.
394 62 : if (dimrho==1) then
395 61 : read(unqps,*)n1,n2,n3
396 :
397 244 : if (all(ngfftf(1:3)== [n1, n2, n3]) ) then
398 61 : read(unqps,*)rhor_out(:,:)
399 : else
400 : write(msg,'(2a,a,5(i3,a),i3)')&
401 0 : 'FFT meshes differ. Performing Fourier interpolation. ',ch10,&
402 0 : 'Found: ',n1,' x',n2,' x',n3,'; Expected: ',ngfftf(1),' x',ngfftf(2),' x',ngfftf(3)
403 0 : ABI_COMMENT(msg)
404 :
405 0 : ABI_MALLOC(rhor_tmp,(n1*n2*n3,nspden))
406 0 : read(unqps,*)rhor_tmp(:,:)
407 :
408 : if (use_FFT_interpolation) then
409 0 : ngfft_found(1:3)=(/n1,n2,n3/)
410 0 : ngfft_found(4)=2*(ngfft_found(1)/2)+1 ! 4:18 are not used, anyway!
411 0 : ngfft_found(5)=2*(ngfft_found(2)/2)+1
412 0 : ngfft_found(6)=ngfft_found(3)
413 0 : ngfft_found(7:18)=ngfftf(7:18)
414 0 : nfft_found=PRODUCT(ngfft_found(1:3)) !no FFT para
415 :
416 0 : cplex_fft =1 ! Real quantities.
417 0 : optin =0 ! Input is taken from rhor.
418 0 : optout =0 ! Output is only in real space.
419 0 : call MPI_enreg%distribfft%free()
420 0 : call MPI_enreg%distribfft%init('c',MPI_enreg%nproc_fft,ngfftf(2),ngfftf(3))
421 0 : call MPI_enreg%distribfft%init('f',MPI_enreg%nproc_fft,ngfft_found(2),ngfft_found(3))
422 :
423 : call fourier_interpol(cplex_fft,nspden,optin,optout,nfft_found,ngfft_found,nfftot,ngfftf,&
424 0 : MPI_enreg,rhor_tmp,rhor_out,rhogdum,rhogdum)
425 :
426 : else
427 : ! Linear interpolation.
428 : do ispden=1,nspden
429 : do ir3=0,ngfftf(3)-1
430 : rr(3)=DBLE(ir3)/n3
431 : do ir2=0,ngfftf(2)-1
432 : rr(2)=DBLE(ir2)/n2
433 : do ir1=0,ngfftf(1)-1
434 : rr(1)=DBLE(ir1)/n1
435 : ifft = 1 +ir1 +ir2*ngfftf(1) +ir3*ngfftf(1)*ngfftf(2)
436 : rhor_out(ifft,ispden) = interpol3d_0d(rr,n1,n2,n3,rhor_tmp(:,ispden))
437 : end do
438 : end do
439 : end do
440 : end do
441 : end if
442 :
443 0 : ABI_FREE(rhor_tmp)
444 : end if
445 :
446 : ! Test the normalization of the QPS density.
447 : ! There might be errors due to the interpolation or the truncation of the G basis set
448 : ! Density will be renormalized in the caller since for PAW we still have to add the onsite contribution.
449 61 : if (usepaw==0) then
450 1806896 : nelect_qps=SUM(rhor_out(:,1))*ucvol/nfftot; ratio=BSt%nelect/nelect_qps
451 : write(msg,'(3(a,f9.4))')&
452 61 : ' Number of electrons calculated using the QPS density = ',nelect_qps,' Expected = ',BSt%nelect,' ratio = ',ratio
453 61 : call wrtout(std_out, msg)
454 : !!rhor_out(:,:)=ratio*rhor_out(:,:)
455 : end if
456 :
457 61 : if (usepaw==1) then
458 : ! Write QP_rhoij for on-site density mixing.
459 0 : read(unqps,*,iostat=ios)natomR,ntypatR
460 0 : if (ios/=0) then
461 0 : msg="Old version of QPS file found. DO NOT USE rhoqpmix for this run."
462 0 : ABI_WARNING(msg)
463 0 : call wrtout(ab_out,msg)
464 : ! Init dummy rhoij just to avoid problems in sigma when rhoij is freed.
465 0 : call pawrhoij_inquire_dim(nspden_rhoij=nspdenR, nspden=nspden)
466 0 : call pawrhoij_alloc(Pawrhoij,1,nspdenR,BSt%nspinor,BSt%nsppol,Cryst%typat,pawtab=Pawtab)
467 0 : close(unqps)
468 0 : RETURN
469 : end if
470 :
471 0 : ABI_CHECK(natomR == Cryst%natom, "mismatch in natom")
472 0 : ABI_CHECK(ntypatR == Cryst%ntypat,"mismatch in ntypat")
473 0 : ABI_MALLOC(nlmn_type, (ntypatR))
474 0 : ABI_MALLOC(typatR, (ntypatR))
475 :
476 0 : read(unqps,*)(typatR(iatom), iatom=1,natomR)
477 0 : ABI_CHECK(ALL(Cryst%typat==typatR), "mismatch in typat")
478 :
479 0 : read(unqps,*)(nlmn_type(itypat), itypat=1,ntypatR)
480 0 : do itypat =1,Cryst%ntypat
481 0 : if (nlmn_type(itypat)/=Pawtab(itypat)%lmn_size) then
482 0 : ABI_ERROR("mismatch in nlmn_type, check QPS file")
483 : end if
484 : end do
485 :
486 0 : read(unqps,*) nsppolR,nspdenR
487 0 : ABI_CHECK(nsppolR==BSt%nsppol,"mismatch in nsppol")
488 0 : ABI_CHECK(nspdenR==nspden ,"mismatch in nspden")
489 :
490 : call pawrhoij_io(pawrhoij,unqps,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,&
491 0 : HDR_LATEST_HEADFORM,"Read",form="formatted")
492 : !% call pawrhoij_io(pawrhoij,std_out,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,HDR_LATEST_HEADFORM,"Echo")
493 :
494 0 : ABI_FREE(nlmn_type)
495 0 : ABI_FREE(typatR)
496 : end if ! usepaw
497 :
498 : end if !dimrho=1
499 :
500 62 : close(unqps)
501 :
502 : else
503 0 : ABI_ERROR("netdf format not implemented")
504 : end if
505 :
506 : DBG_EXIT("COLL")
507 :
508 126 : end subroutine rdqps
509 : !!***
510 :
511 : !----------------------------------------------------------------------
512 :
513 : !!****f* m_qparticles/show_QP
514 : !! NAME
515 : !! show_QP
516 : !!
517 : !! FUNCTION
518 : !! Print in a nice format (?) the expansion coefficients of the quasiparticle
519 : !! amplitudes in terms of the KS eigenvectors
520 : !!
521 : !! INPUTS
522 : !! Bst<ebands_t>=Description of the band structure.
523 : !! %nsppol=1 for unpolarized, 2 for spin-polarized.
524 : !! %mband=Max number of bands (in GW doesn"t depend on k an spin)
525 : !! %nkpt=number of irreducible k-points.
526 : !! %eig(mband,nkpt,nsppol)= QP energies for each k-point, band and spin.
527 : !! m_ks_to_qp(nbnds,nbnds,nkibz,nsppol)=matrix giving the decomposition of the QP
528 : !! amplitued in the mainfold generated by the KS wavefunctions
529 : !! (i.e $ m_ks_to_qp(ib,jb,k,s) := \langle \psi_{ib,k,s}^{KS}| \psi_{jb,k,s}^{QP}\rangle $
530 : !! fromb,tob=initial and final band index for QP, only states in this range are printed
531 : !! prtvol=Verbosity level (not used)
532 : !! unit=Unit number of the output file
533 : !! tolmat[Optional]=Only components whose coefficient has modulus larger than tolmat are shown (default is 0.01)
534 : !!
535 : !! OUTPUT
536 : !! Only printing
537 : !!
538 : !! NOTES
539 : !! Only master node should call this routine.
540 : !!
541 : !! SOURCE
542 :
543 0 : subroutine show_QP(Bst,m_ks_to_qp,fromb,tob,unit,prtvol,tolmat,kmask)
544 :
545 : !Arguments ------------------------------------
546 : !scalars
547 : integer,optional,intent(in) :: fromb,tob
548 : integer,optional,intent(in) :: prtvol,unit
549 : real(dp),optional,intent(in) :: tolmat
550 : type(ebands_t),intent(in) :: Bst
551 : !arrays
552 : logical,optional,intent(in) :: kmask(Bst%nkpt)
553 : complex(dp),intent(in) :: m_ks_to_qp(Bst%mband,Bst%mband,Bst%nkpt,Bst%nsppol)
554 :
555 : !Local variables-------------------------------
556 : !scalars
557 : integer,parameter :: NBRA=5
558 : logical,parameter :: use_rhophi=.True.
559 : integer :: ib_start,ib_stop,my_prtvol,counter,ib_KS,ib_QP,ikibz,isp,nspace,my_unt,nband_k
560 : real(dp) :: my_tolmat,rho,phi
561 : character(len=10) :: bks,bqp,k_tag,spin_tag
562 : character(len=500) :: KS_row,KS_ket,tmpstr,QP_ket
563 : !arrays
564 : real(dp) :: cx(2)
565 : ! *********************************************************************
566 :
567 28 : my_unt =std_out ; if (PRESENT(unit )) my_unt =unit
568 28 : my_prtvol=0 ; if (PRESENT(prtvol)) my_prtvol=prtvol
569 28 : ib_start =1 ; if (PRESENT(fromb )) ib_start =fromb
570 28 : ib_stop =Bst%mband; if (PRESENT(tob )) ib_stop =tob
571 28 : my_tolmat=0.001 ; if (PRESENT(tolmat)) my_tolmat=ABS(tolmat)
572 :
573 : ! I suppose nband_k is constant thus the check is done here.
574 28 : if (ib_start<=0 ) ib_start=1
575 28 : if (ib_start>Bst%mband) ib_start=Bst%mband
576 28 : if (ib_stop<=0 ) ib_stop=1
577 28 : if (ib_stop>Bst%mband ) ib_stop=Bst%mband
578 :
579 : ! Have to follow rules 7.f.
580 : write(my_unt,'(/,a,/,a,/,a,f6.3,a,/,a)')&
581 28 : ' '//REPEAT('*',76),&
582 28 : & ' ***** QP amplitudes expressed as linear combination of KS eigenstates. *****',&
583 28 : & ' ***** Only KS components whose modulus is larger than ',my_tolmat,' are shown ***** ',&
584 56 : & ' '//REPEAT('*',76)
585 : if (use_rhophi) then
586 28 : write(my_unt,"(a)")"Complex coefficients given in (rho, phi) polar representation."
587 : else
588 : write(my_unt,"(a)")"Complex coefficients given in (Re, Im) representation."
589 : end if
590 :
591 28 : if (PRESENT(kmask)) then
592 0 : if (.not.ALL(kmask)) write(my_unt,'(/,a,i3,a)')' Only ',COUNT(kmask),' k-points are reported '
593 : end if
594 :
595 56 : do isp=1,Bst%nsppol
596 28 : call int2char10(isp,spin_tag)
597 28 : write(my_unt,'(/,a,i2,a,/)')' >>>>> Begin block for spin ',isp,' <<<<< '
598 :
599 232 : do ikibz=1,Bst%nkpt
600 176 : if (PRESENT(kmask)) then
601 0 : if (.not.kmask(ikibz)) CYCLE
602 : end if
603 176 : call int2char10(ikibz,k_tag)
604 176 : nband_k=Bst%nband(ikibz+(isp-1)*Bst%nkpt)
605 176 : write(my_unt,'(a,i4,a,3es16.8,a,f6.3,/)')' k-point: ',ikibz,') ',Bst%kptns(:,ikibz),'; wtk= ',Bst%wtk(ikibz)
606 :
607 2288 : do ib_QP=ib_start,ib_stop
608 2084 : call int2char10(ib_QP,bqp)
609 2084 : QP_ket=' |QP: b='//TRIM(bqp)//'; s='//TRIM(spin_tag)//'> = '
610 2084 : write(my_unt,'(a)')TRIM(QP_ket)
611 2084 : nspace=LEN(TRIM(QP_ket))
612 :
613 2084 : counter=0 ; KS_row=REPEAT('',nspace+2)
614 30076 : do ib_KS=1,Bst%mband
615 27992 : if (ABS(m_ks_to_qp(ib_KS,ib_QP,ikibz,isp))<my_tolmat) CYCLE
616 8743 : counter=counter+1
617 8743 : call int2char10(ib_KS,bks)
618 8743 : write(tmpstr,'(3a)')' |',TRIM(bks),'>'
619 :
620 : if (use_rhophi) then
621 : ! coefficient as (rho, phi)
622 8743 : cx(1) = real(m_ks_to_qp(ib_KS,ib_QP,ikibz,isp))
623 8743 : cx(2) = aimag(m_ks_to_qp(ib_KS,ib_QP,ikibz,isp))
624 8743 : call rhophi(cx, phi, rho)
625 8743 : write(KS_ket,'(1x,2f7.3,a,1x)')rho, phi, TRIM(tmpstr)
626 : else
627 : ! coefficient as (Re, Im)
628 : write(KS_ket,'(1x,2f7.3,a,1x)')m_ks_to_qp(ib_KS,ib_QP,ikibz,isp),TRIM(tmpstr)
629 : end if
630 8743 : KS_row=TRIM(KS_row)//TRIM(KS_ket)
631 19570 : if (MOD(counter,NBRA)==0) then ! nbra KS kets per row
632 821 : write(my_unt,'(a)')TRIM(KS_row)
633 821 : KS_row=REPEAT('',nspace+2)
634 : end if
635 : end do
636 :
637 2084 : if (MOD(counter,NBRA)/=0) write(my_unt,'(a)')TRIM(KS_row) ! Last row, if any
638 2260 : write(my_unt,'(a)')''
639 : end do !ib_QP
640 :
641 : end do !ikibz
642 : end do !isp
643 :
644 28 : write(my_unt,'(a,/)')' '//REPEAT('*',76)
645 :
646 28 : end subroutine show_QP
647 : !!***
648 :
649 : !----------------------------------------------------------------------
650 :
651 : !!****f* m_qparticles/rdgw
652 : !! NAME
653 : !! rdgw
654 : !!
655 : !! FUNCTION
656 : !! This subroutine reads the GW corrections from a _GW file.
657 : !!
658 : !! INPUTS
659 : !! [extrapolate]= if .TRUE., the routine extrapolates the
660 : !! GW corrections for the states that have not been explicitly evaluated (default).
661 : !! If .FALSE., only the GW states that have been calculated will be used to replace
662 : !! the input eigenvalues stored in Bst%eig
663 : !! Bst<ebands_t>=type describing the Band structure.
664 : !! %nbnds=number of bands.
665 : !! %nkpt=number of irred k-points.
666 : !! %nsppol=number of spin
667 : !! %kptns(3,nkpt)=irreducible k-points
668 : !!
669 : !! SIDE EFFECTS
670 : !! Bst%eig(%mband,%nkpt,%nsppol)=Overwritten with GW energies according to extrapolate flag.
671 : !!
672 : !! OUTPUT
673 : !! igwene(Bst%mband,Bst%nkpt,Bst%nsppol)= The imaginary part of the QP energies.
674 : !!
675 : !! SOURCE
676 :
677 0 : subroutine rdgw(Bst,fname,igwene,extrapolate)
678 :
679 : !Arguments ------------------------------------
680 : !scalars
681 : character(len=*),intent(in) :: fname
682 : logical,optional,intent(in) :: extrapolate
683 : type(ebands_t),intent(inout) :: Bst
684 : !arrays
685 : real(dp),intent(out) :: igwene(Bst%mband,Bst%nkpt,Bst%nsppol)
686 :
687 : !Local variables ------------------------------
688 : !scalars
689 : integer :: ib,ibr,ik,ikibz,ikr,is,nn,nbandR,nkibzR,nsppolR,unt,nbv
690 : real(dp) :: alpha,beta,degw,egw_r,egw_i,smrt
691 : logical :: do_extrapolate
692 : character(len=500) :: msg
693 : !arrays
694 0 : integer,allocatable :: vbik(:,:),seen(:)
695 : real(dp) :: kread(3)
696 0 : real(dp),allocatable :: gwcorr(:,:,:)
697 : !************************************************************************
698 :
699 0 : call wrtout(std_out,'Reading GW corrections from file: '//TRIM(fname))
700 0 : ABI_CHECK(ALL(Bst%nband==Bst%mband),"nband must be constant")
701 :
702 0 : if (open_file(fname,msg,newunit=unt,status='old') /=0) then
703 0 : ABI_ERROR(msg)
704 : end if
705 :
706 0 : read(unt,*)nkibzR,nsppolR
707 :
708 0 : ABI_CHECK(nsppolR==Bst%nsppol,"mismatch in nsppol")
709 0 : if (nkibzR/=Bst%nkpt) then
710 : write(msg,'(a,i4,a,i4,2a)')&
711 0 : 'Found less k-points than that required ',nkibzR,'/',Bst%nkpt,ch10,&
712 0 : 'Some k-points will be skipped. Continuing anyway '
713 0 : ABI_WARNING(msg)
714 : end if
715 :
716 0 : ABI_MALLOC(gwcorr,(Bst%mband,Bst%nkpt,Bst%nsppol))
717 0 : ABI_MALLOC(seen,(Bst%nkpt))
718 0 : gwcorr=zero
719 0 : igwene=zero
720 :
721 0 : do is=1,Bst%nsppol
722 0 : seen=0
723 :
724 0 : do ikr=1,nkibzR
725 0 : read(unt,*)kread(:)
726 0 : read(unt,*)nbandR
727 0 : ikibz=0
728 0 : do ik=1,Bst%nkpt
729 0 : if (ALL(ABS(kread(:)-Bst%kptns(:,ik))<0.0001)) then
730 0 : ikibz=ik
731 0 : seen(ik) = seen(ik) + 1
732 : end if
733 : end do
734 0 : do ib=1,nbandR
735 0 : read(unt,*)ibr,egw_r,degw,egw_i
736 0 : if (ibr<=Bst%mband .and. ikibz/=0) then
737 0 : gwcorr(ibr,ikibz,is)=degw/Ha_eV
738 0 : igwene(ibr,ikibz,is)=egw_i/Ha_eV
739 : end if
740 : end do
741 : end do
742 :
743 0 : if (ANY(seen/=1)) then
744 0 : do ik=1,Bst%nkpt
745 0 : if (seen(ik)/=1) then
746 0 : write(msg,'(a,3f8.3,a)')" k-point: ",Bst%kptns(:,ik)," not found in the GW file!"
747 0 : ABI_WARNING(msg)
748 : end if
749 : end do
750 : end if
751 :
752 : end do
753 :
754 0 : ABI_FREE(seen)
755 0 : close(unt)
756 :
757 0 : do_extrapolate=.TRUE.; if (PRESENT(extrapolate)) do_extrapolate=extrapolate
758 :
759 0 : if (.not. do_extrapolate) then ! Only the bands calculated are updated.
760 0 : Bst%eig = Bst%eig + gwcorr
761 :
762 : else
763 :
764 0 : if (ANY(ABS(igwene)>tol6)) then
765 0 : write(msg,'(4a)')ch10,&
766 0 : "The GW file contains QP energies with non-zero imaginary part",ch10,&
767 0 : "Extrapolation not coded, change the source! "
768 0 : ABI_ERROR(msg)
769 : end if
770 :
771 0 : ABI_MALLOC(vbik,(BSt%nkpt,BSt%nsppol))
772 0 : vbik(:,:) = BSt%get_valence_idx()
773 :
774 0 : do is=1,Bst%nsppol
775 0 : do ik=1,Bst%nkpt
776 :
777 0 : nbv=vbik(ik,is) ! Index of the (valence band| Fermi band) for each spin
778 0 : nn=Bst%mband-nbv
779 :
780 0 : do ib=nbv+1,Bst%mband
781 0 : if ( ABS(gwcorr(ib,ik,is)) < tol16) then
782 0 : nn=ib-1-nbv
783 0 : if (nn>1) then
784 0 : call wrtout(std_out,"Linear extrapolating (conduction) GW corrections beyond the read values")
785 0 : smrt=linfit(nn,Bst%eig(nbv+1:nbv+nn,ik,is),gwcorr(nbv+1:nbv+nn,ik,is),alpha,beta)
786 : else
787 0 : call wrtout(std_out,"Assuming constant (conduction) GW corrections beyond the read values")
788 0 : alpha=zero
789 0 : beta =gwcorr(nbv+nn,ik,is)
790 : end if
791 : EXIT !ib loop
792 : end if
793 : end do !ib
794 :
795 0 : do ib=nbv+nn+1,Bst%mband
796 0 : gwcorr(ib,ik,is)= alpha*Bst%eig(ib,ik,is) + beta
797 : end do
798 :
799 0 : nn=nbv
800 0 : do ib=nbv,1,-1
801 0 : if ( ABS(gwcorr(ib,ik,is)) < tol16) then
802 0 : nn=nbv-ib
803 0 : if (nn>1) then
804 0 : call wrtout(std_out,"Linear extrapolating (valence) GW corrections beyond the read values")
805 0 : smrt=linfit(nn,Bst%eig(nbv-nn+1:nbv,ik,is),gwcorr(nbv-nn+1:nbv,ik,is),alpha,beta)
806 : else
807 0 : call wrtout(std_out,"Assuming constant (valence) GW corrections beyond the read values")
808 0 : alpha=zero
809 0 : beta =gwcorr(nbv,ik,is)
810 : end if
811 : EXIT !ib
812 : end if
813 : end do !ib
814 :
815 0 : do ib=1,nbv-nn
816 0 : gwcorr(ib,ik,is)=alpha*Bst%eig(ib,ik,is) + beta
817 : end do
818 :
819 : end do !ik
820 : end do !is
821 :
822 0 : call wrtout(std_out,' k s GW corrections [eV] ')
823 0 : do is=1,Bst%nsppol
824 0 : do ik=1,Bst%nkpt
825 0 : write(msg,'(i3,1x,i3,10f7.2/50(10x,10f7.2/))')ik,is,(Ha_eV*gwcorr(ib,ik,is),ib=1,Bst%mband)
826 0 : call wrtout(std_out,msg)
827 : end do
828 : end do
829 0 : Bst%eig = Bst%eig + gwcorr
830 0 : ABI_FREE(vbik)
831 : end if
832 :
833 0 : call wrtout(std_out,' k s GW eigenvalues [eV]')
834 0 : do is=1,Bst%nsppol
835 0 : do ik=1,Bst%nkpt
836 0 : write(std_out,'(2(i3,1x),7x,10f7.2/50(15x,10f7.2/))')ik,is,(Ha_eV*Bst%eig(ib,ik,is),ib=1,Bst%mband)
837 : end do
838 : end do
839 :
840 0 : ABI_FREE(gwcorr)
841 :
842 0 : end subroutine rdgw
843 : !!***
844 :
845 : !----------------------------------------------------------------------
846 :
847 : !!****f* m_qparticles/updt_m_ks_to_qp
848 : !! NAME
849 : !! updt_m_ks_to_qp
850 : !!
851 : !! FUNCTION
852 : !! Updates the matrix containing the unitary transformation from the lda states
853 : !! to the quasiparticle states.
854 : !!
855 : !! INPUTS
856 : !! Sigp<sigparams_t>=Parameters characterizing the self-energy calculation.
857 : !! %nsppol=1 for unpolarized, 2 for spin-polarized
858 : !! %nbnds=number of bands used for sigma
859 : !! Sr<sigma_t>=Structure containing the results of the sigma run.
860 : !! %en_qp_diago(nbnds,nibz,nsppol)= NEW quasi-particle energies
861 : !! %eigvec_qp(nbnds,nbnds,nibz,nsppol)= NEW QP amplitudes in the KS basis set
862 : !! obtained by diagonalizing H0 + Herm(Sigma).
863 : !! Kmesh<kmesh_t>=information on the k-point sampling.
864 : !! %nibz=number of irreducible k-points
865 : !! %ibz(3,kibz)=reduced coordinates of the irreducible k-points
866 : !! nscf=Number of self consistent cycles performed
867 : !!
868 : !! OUTPUT
869 : !! (see side effects)
870 : !!
871 : !! SIDE EFFECTS
872 : !! m_ks_to_qp(nbnds,nbnds,nibz,nsppol)= overwritten with the new QP amplitudes
873 : !! in terms of KS wavefunctions
874 : !!
875 : !! NOTES
876 : !! Only master node should call this routine.
877 : !!
878 : !! SOURCE
879 :
880 71 : subroutine updt_m_ks_to_qp(Sigp,Kmesh,nscf,Sr,m_ks_to_qp)
881 :
882 : !Arguments ------------------------------------
883 : !scalars
884 : integer,intent(in) :: nscf
885 : type(kmesh_t),intent(in) :: Kmesh
886 : type(sigparams_t),intent(in) :: Sigp
887 : type(sigma_t),intent(in) :: Sr
888 : !arrays
889 : complex(dp),intent(inout) :: m_ks_to_qp(Sigp%nbnds,Sigp%nbnds,Kmesh%nibz,Sigp%nsppol)
890 :
891 : !Local variables-------------------------------
892 : !scalars
893 : integer :: ik,is
894 : !arrays
895 71 : complex(dp),allocatable :: mtmp(:,:)
896 : ! *************************************************************************
897 :
898 71 : if (nscf >= 0) then
899 : ! Calculate the new m_ks_to_qp
900 71 : ABI_CHECK(Sr%needs_eigvec_qp, "needs_eigvec should be true")
901 284 : ABI_MALLOC(mtmp, (Sigp%nbnds,Sigp%nbnds))
902 144 : do is=1,Sigp%nsppol
903 567 : do ik=1,Kmesh%nibz
904 71003 : mtmp(:,:)=m_ks_to_qp(:,:,ik,is)
905 2086190 : m_ks_to_qp(:,:,ik,is)=MATMUL(mtmp(:,:),Sr%eigvec_qp(:,:,ik,is))
906 : end do
907 : end do
908 71 : ABI_FREE(mtmp)
909 : end if
910 :
911 71 : end subroutine updt_m_ks_to_qp
912 :
913 : !----------------------------------------------------------------------
914 :
915 : END MODULE m_qparticles
916 : !!***
|
