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