Line data Source code
1 : !!****m* ABINIT/m_gwrdm
2 : !! NAME
3 : !! m_gwrdm
4 : !!
5 : !! FUNCTION
6 : !! Compute density matrix correction Galitskii-Migdal Ecorr, G = Go + Go Sigma Go (imaginary freqs. are used in Sigma_c)
7 : !! and associated quantities (natural orbitals, matrix elements, etc.).
8 : !!
9 : !! SOURCE
10 :
11 : #if defined HAVE_CONFIG_H
12 : #include "config.h"
13 : #endif
14 :
15 : #include "abi_common.h"
16 :
17 : module m_gwrdm
18 :
19 : use defs_basis
20 : use m_gwdefs
21 : use m_abicore
22 : use m_xmpi
23 : use m_errors
24 : use m_hide_blas
25 : use m_time
26 : use m_wfd
27 : use m_hdr
28 : use m_dtset
29 :
30 : use m_fstrings, only : sjoin, itoa
31 : use m_melemts, only : melements_t
32 : use m_bz_mesh, only : kmesh_t
33 : use m_ebands, only : ebands_t
34 : use m_sigma, only : sigma_t
35 : use m_xctk, only : xcden
36 : use m_gaussian_quadrature, only: cgqf
37 :
38 : implicit none
39 :
40 : private :: no2ks,ks2no,printrdm_k,rotate_ks_no
41 : !!***
42 :
43 : public :: quadrature_sigma_cw,calc_Ec_GM_k,calc_rdmx,calc_rdmc,natoccs,update_hdr_bst,print_tot_occ,change_matrix
44 : public :: print_chkprdm,print_total_energy,print_band_energies,get_chkprdm
45 : !!***
46 :
47 : contains
48 : !!***
49 :
50 : !!****f* ABINIT/quadrature_sigma_cw
51 : !! NAME
52 : !! quadrature_sigma_cw
53 : !!
54 : !! FUNCTION
55 : !! Quadrature frequencies used for Sigma_c(iw) integration
56 : !!
57 : !! INPUTS
58 : !! Sigp<sigparams_t>=Parameters governing the self-energy calculation.
59 : !! Sr=sigma_t (see the definition of this structured datatype)
60 : !! weights=real quadrature weights.
61 : !!
62 : !! OUTPUT
63 : !! Update Sigp and Sr imaginary frequencies with iw, and weights with the quadrature weights
64 : !!
65 : !! SOURCE
66 :
67 5 : subroutine quadrature_sigma_cw(Sigp,Sr,weights)
68 : !Arguments ------------------------------------
69 : !scalars
70 : type(sigparams_t),intent(inout) :: Sigp
71 : type(sigma_t),intent(inout) :: Sr
72 : !arrays
73 : real(dp),intent(inout) :: weights(:)
74 :
75 : !Local variables ------------------------------
76 : !scalars
77 : integer :: ifreqs,order_int,gaussian_kind,units(2)
78 : real(dp) :: gwalpha,gwbeta,wmin,wmax
79 : character(len=500) :: msg
80 : !arrays
81 5 : real(dp),allocatable :: freqs(:)
82 : !************************************************************************
83 :
84 15 : units = [std_out, ab_out]
85 :
86 5 : order_int=Sigp%nomegasi
87 5 : write(msg,'(a45,i9)')' number of imaginary frequencies for Sigma_c ',order_int
88 5 : call wrtout(units, msg)
89 5 : write(msg,'(a1)')' '
90 5 : call wrtout(units, msg)
91 5 : order_int=Sigp%nomegasi
92 15 : ABI_MALLOC(freqs,(order_int))
93 5 : gaussian_kind=1
94 5 : gwalpha=zero
95 5 : gwbeta=zero
96 5 : wmin=zero
97 5 : wmax=one
98 5 : call cgqf(order_int,gaussian_kind,gwalpha,gwbeta,wmin,wmax,freqs,weights)
99 : ! From 0 to 1 -> 0 to infinity
100 165 : weights(:)=weights(:)/(one-freqs(:))**two
101 165 : freqs(:)=freqs(:)/(one-freqs(:))
102 : ! Form complex frequencies from 0 to iInf and print them in the log file
103 5 : write(msg,'(a52)')' Re(iw) Im(iw) Weight '
104 5 : call wrtout(std_out,msg)
105 5 : write(msg,'(a52)')' -------- -------- -------- '
106 5 : call wrtout(std_out,msg)
107 165 : do ifreqs=1,order_int
108 160 : Sigp%omegasi(ifreqs)=cmplx(zero,freqs(ifreqs))
109 160 : Sr%omega_i(ifreqs)=Sigp%omegasi(ifreqs)
110 160 : write(msg,'(3f17.5)') Sr%omega_i(ifreqs),weights(ifreqs)
111 165 : call wrtout(std_out,msg)
112 : enddo
113 5 : ABI_FREE(freqs)
114 :
115 5 : end subroutine quadrature_sigma_cw
116 : !!***
117 :
118 : !!****f* ABINIT/Calc_Ec_GM_k
119 : !! NAME
120 : !! calc_Ec_GM_k
121 : !!
122 : !! FUNCTION
123 : !! Calculate Galitskii-Migdal corr. energy integrated in the Imaginary axis Ec = 1/pi \sum_i \int Gii(iv)*Sigma_c,ii(iv) + cc. dv
124 : !!
125 : !! INPUTS
126 : !! ib1=min band for given k
127 : !! ib2=max band for given k.
128 : !! ik_ibz= the label of k-point in the IBZ whose Galitskii-Migdal contribution is accounted.
129 : !! weights=array containing the weights used in the quadrature.
130 : !! sigcme_k=array containing Sigma(iw) as Sigma(iw,ib1:ib2,ib1:ib2,nspin)
131 : !! rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details).
132 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
133 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
134 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
135 : !! Sr=sigma_t (see the definition of this structured datatype)
136 : !!
137 : !! OUTPUT
138 : !! Compute the Galitskii-Migdal corr energy contribution of this k-point:
139 : !! Ec ^k = 1/(4*pi) * fact_spin * int _{ -Inf }^{ +Inf } dv Sigma_c ^k (iv) * G0(iv)
140 : !! = 1/(4*pi) * fact_spin * int _{ 0 }^{ +Inf } dv 2 * Re{ Sigma_c ^k (iv) * G0(iv) }
141 : !!
142 : !! SOURCE
143 :
144 0 : function calc_Ec_GM_k(ib1,ib2,ik_ibz,Sr,weights,sigcme_k,ebands) result(Ec_GM_k)
145 :
146 : !Arguments ------------------------------------
147 : !scalars
148 : real(dp) :: Ec_GM_k
149 : integer,intent(in) :: ib1,ib2,ik_ibz
150 : type(ebands_t),target,intent(in) :: ebands
151 : type(sigma_t),intent(in) :: Sr
152 : !arrays
153 : real(dp),intent(in) :: weights(:)
154 : complex(dp),intent(in) :: sigcme_k(:,:,:,:)
155 : !Local variables ------------------------------
156 : !scalars
157 : integer :: ibdm!,unitt
158 : real(dp) :: ec_integrated,spin_fact,fact
159 : character(len=500) :: msg
160 : !arrays
161 : !************************************************************************
162 :
163 0 : ec_integrated=zero
164 0 : spin_fact=two
165 0 : fact=spin_fact*(one/(two_pi*two))
166 :
167 0 : if (ib1/=1) then
168 : msg="Unable to compute the Galitskii-Migdal correlation energy because the first band was " // &
169 0 : & "not included in bdgw interval. Restart the calculation starting bdgw from 1."
170 0 : ABI_WARNING(msg)
171 : else
172 : ! WARNING: Sigma_c(iv) produced from a previous integration at the screening stage, is numerically not much stable and introduces bumps.
173 : ! Unfortunately, the Green's function times Sigma_c(iv) does not decay fast enough with iv to overcome the bumps. These bumps are
174 : ! not pronouced for the linearized density matrix update, as two Green's functions are multiplied making the decay much faster with iv.
175 : ! If a better way to produce more stable Sigma_c(iv) values is found, this subroutine can be use to evaluate GM Ecorr in the future. TODO
176 0 : do ibdm=1,ib2
177 : ! Sigma_pp(iv)/[(iv - e_ibdm,k)] + [Sigma_pp(iv)/[(iv - e_ibdm,k)]]^* = 2 Re [Sigma_pp(iv)/(iv - e_ibdm,k)]
178 0 : ec_integrated=ec_integrated+two*real( sum(weights(:)*sigcme_k(:,ibdm,ibdm,1)/(Sr%omega_i(:)-ebands%eig(ibdm,ik_ibz,1)) ) )
179 : end do
180 : endif
181 :
182 0 : Ec_GM_k=fact*ec_integrated
183 :
184 0 : end function calc_Ec_GM_k
185 : !!***
186 :
187 : !!****f* ABINIT/calc_rdmx
188 : !! NAME
189 : !! calc_rdmx
190 : !!
191 : !! FUNCTION
192 : !! Calculate density matrix corrections for G = Go + Go (Sigma_x - alpha*Sigma_x - Vxc) Go
193 : !!
194 : !! INPUTS
195 : !! ib1=min band for given k
196 : !! ib2=max band for given k.
197 : !! ik_ibz= the label of k-point in the IBZ.
198 : !! rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details).
199 : !! pot=Self-energy-Potential difference, matrix size (i,j), where i and j belong to k.
200 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
201 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
202 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
203 : !!
204 : !! OUTPUT
205 : !! Updated rdm_k matrix array with Go (Sigma_x - alpha*Sigma_x - Vxc) Go
206 : !!
207 : !! SOURCE
208 :
209 27 : subroutine calc_rdmx(ib1,ib2,ik_ibz,pot,rdm_k,ebands)
210 :
211 : !Arguments ------------------------------------
212 : !scalars
213 : integer,intent(in) :: ib1,ib2,ik_ibz
214 : type(ebands_t),target,intent(in) :: ebands
215 : !arrays
216 : complex(dp),intent(in) :: pot(:,:)
217 : complex(dp),intent(inout) :: rdm_k(:,:)
218 :
219 : !Local variables ------------------------------
220 : !scalars
221 : character(len=500) :: msg
222 : integer :: ib1dm,ib2dm, units(2)
223 : real(dp) :: spin_fact,tol8
224 : !************************************************************************
225 :
226 27 : tol8=1.0e-8
227 27 : spin_fact=two
228 81 : units = [std_out, ab_out]
229 :
230 27 : write(msg,'(a58,3f10.5)')' Computing the 1-RDM correction for Sx-Vxc and k-point: ',ebands%kptns(:,ik_ibz)
231 27 : call wrtout(units, msg)
232 27 : write(msg,'(a11,i5,a8,i5)')'from band ',ib1,' to band',ib2
233 27 : call wrtout(units, msg)
234 :
235 1709 : rdm_k(:,:)=czero
236 196 : do ib1dm=ib1,ib2-1
237 841 : do ib2dm=ib1dm+1,ib2
238 814 : if ((ebands%occ(ib1dm,ik_ibz,1)>tol8) .and. (ebands%occ(ib2dm,ik_ibz,1)<tol8)) then
239 : rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1))=spin_fact&
240 236 : &*pot(1+(ib1dm-ib1),1+(ib2dm-ib1))/(ebands%eig(ib1dm,ik_ibz,1)-ebands%eig(ib2dm,ik_ibz,1)+tol8)
241 : ! Dji = Dij^*
242 236 : rdm_k(1+(ib2dm-ib1),1+(ib1dm-ib1))=conjg(rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1)))
243 : end if
244 : end do
245 : end do
246 :
247 27 : end subroutine calc_rdmx
248 : !!***
249 :
250 : !!****f* ABINIT/calc_rdmc
251 : !! NAME
252 : !! calc_rdmc
253 : !!
254 : !! FUNCTION
255 : !! Calculate density matrix corrections for G = Go + int Go(iw) Sigma_c(iw) Go(iw) dw
256 : !!
257 : !! INPUTS
258 : !! ib1=min band for given k
259 : !! ib2=max band for given k.
260 : !! ik_ibz= the label of k-point in the IBZ.
261 : !! omega_i=Frequencies along the imaginary axis.
262 : !! weights=array containing the weights used in the quadrature.
263 : !! sigcme_k=array containing Sigma(iw) as Sigma(iw,ib1:ib2,ib1:ib2,nspin)
264 : !! rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details).
265 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
266 : !!
267 : !! OUTPUT
268 : !! Updated rdm_k matrix array with int Go(iw) Sigma_c(iw) Go(iw) dw
269 : !!
270 : !! SOURCE
271 :
272 21 : subroutine calc_rdmc(ib1,ib2,ik_ibz,omega_i,weights,sigcme_k,ebands,rdm_k)
273 :
274 : !Arguments ------------------------------------
275 : !scalars
276 : integer,intent(in) :: ib1,ib2,ik_ibz
277 : type(ebands_t),target,intent(in) :: ebands
278 : complex(dp),intent(in) :: omega_i(:)
279 : !arrays
280 : real(dp),intent(in) :: weights(:)
281 : complex(dp),intent(inout) :: rdm_k(:,:)
282 : complex(dp),intent(in) :: sigcme_k(:,:,:,:)
283 :
284 : !Local variables ------------------------------
285 : !scalars
286 : real(dp) :: spin_fact,fact
287 : integer :: ib1dm, ib2dm, units(2)
288 : character(len=500) :: msg
289 : !************************************************************************
290 :
291 21 : spin_fact = two
292 21 : fact = spin_fact * (one/two_pi)
293 63 : units = [std_out, ab_out]
294 :
295 21 : write(msg,'(a58,3f10.5)')' Computing the 1-RDM correction for Sc(iw) and k-point: ',ebands%kptns(1:,ik_ibz)
296 21 : call wrtout(units , msg)
297 21 : write(msg,'(a11,i5,a8,i5)')'from band ',ib1,' to band',ib2
298 21 : call wrtout(units , msg)
299 :
300 1271 : rdm_k(:,:)=czero
301 169 : do ib1dm=ib1,ib2
302 794 : do ib2dm=ib1dm,ib2
303 : ! Sigma_pq/[(denominator)] + [Sigma_qp/[(denominator)]]^*
304 : rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1))=fact*sum(weights(:)*( sigcme_k(:,1+(ib1dm-ib1),1+(ib2dm-ib1),1)/&
305 : &( (omega_i(:)-ebands%eig(ib1dm,ik_ibz,1))*(omega_i(:)-ebands%eig(ib2dm,ik_ibz,1)) )&
306 : +conjg( sigcme_k(:,1+(ib2dm-ib1),1+(ib1dm-ib1),1)/&
307 19605 : &( (omega_i(:)-ebands%eig(ib1dm,ik_ibz,1))*(omega_i(:)-ebands%eig(ib2dm,ik_ibz,1)) ) ) ) )
308 : ! Dji = Dij^*
309 773 : rdm_k(1+(ib2dm-ib1),1+(ib1dm-ib1))=conjg(rdm_k(1+(ib1dm-ib1),1+(ib2dm-ib1)))
310 : end do
311 : end do
312 :
313 21 : end subroutine calc_rdmc
314 : !!***
315 :
316 : !!****f* ABINIT/natoccs
317 : !! NAME
318 : !! natoccs
319 : !!
320 : !! FUNCTION
321 : !! Calculate natural orbitals and occ. numbers for a given k-point
322 : !!
323 : !! INPUTS
324 : !! ib1=min band for given k
325 : !! ib2=max band for given k.
326 : !! ik_ibz= the label of k-point in the IBZ.
327 : !! iinfo=use Sigma_x or Sigma_c phaser
328 : !! weights=array containing the weights used in the quadrature.
329 : !! nateigv=array containing the natural eigenvectors in columns (nbands,nband,k-point,nspin)
330 : !! rdm_k=density matrix, matrix (i,j), where i and j belong to the k-point k (see m_sigma_driver for more details).
331 : !! occs = array containing the occ numbers for a given k-point occs(nband,k-point).
332 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
333 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
334 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
335 : !! checksij=check the orthonormality of the nat. orbitals
336 : !!
337 : !! OUTPUT
338 : !! Compute the nat. orbitals and occ. numbers from the rdm_k matrix (for exchange and correlations)
339 : !!
340 : !! SOURCE
341 :
342 48 : subroutine natoccs(ib1,ib2,rdm_k,nateigv,occs,ebands,ik_ibz,iinfo,checksij)
343 : !Arguments ------------------------------------
344 : !scalars
345 : integer,intent(in) :: ib1,ib2,ik_ibz,iinfo
346 : integer,intent(in),optional :: checksij
347 : type(ebands_t),target,intent(in) :: ebands
348 : !arrays
349 : real(dp),intent(inout) :: occs(:,:)
350 : complex(dp),intent(inout) :: rdm_k(:,:),nateigv(:,:,:,:)
351 : !Local variables ------------------------------
352 : !scalars
353 : integer:: ndim,ib1dm,ib2dm,ib3dm,lwork,info
354 : logical:: check_Sijmat
355 : character(len=500) :: msg
356 : real(dp) :: toccs_k,tol10
357 : complex(dp) :: Sib1k_ib2k
358 : !arrays
359 : integer :: units(2)
360 48 : real(dp),allocatable :: occs_tmp(:),occs_tmp2(:),rwork(:)
361 48 : complex(dp),allocatable :: work(:),tmp_mat(:,:),eigenvect(:,:)
362 : !************************************************************************
363 :
364 48 : check_Sijmat=.false.; if (present(checksij)) check_Sijmat=.true.
365 144 : units = [std_out, ab_out]
366 48 : tol10=1.0e-10
367 :
368 48 : ndim=ib2-ib1+1
369 48 : lwork=2*ndim-1
370 144 : ABI_MALLOC(occs_tmp,(ndim))
371 96 : ABI_MALLOC(occs_tmp2,(ndim))
372 144 : ABI_MALLOC(work,(lwork))
373 192 : ABI_MALLOC(tmp_mat,(ndim,ndim))
374 144 : ABI_MALLOC(eigenvect,(ndim,ndim))
375 144 : ABI_MALLOC(rwork,(3*ndim-2))
376 :
377 2980 : tmp_mat=zero
378 392 : do ib2dm=1,ndim
379 1858 : do ib1dm=ib2dm,ndim
380 1466 : tmp_mat(ib1dm,ib2dm)=rdm_k(ib1dm,ib2dm)
381 : ! Dji = Dij^*
382 1810 : tmp_mat(ib2dm,ib1dm)=conjg(tmp_mat(ib1dm,ib2dm))
383 : end do
384 : end do
385 :
386 688 : work=zero
387 392 : occs_tmp=zero
388 48 : info=0
389 48 : call zheev('v','u',ndim,tmp_mat,ndim,occs_tmp,work,lwork,rwork,info)
390 48 : ABI_CHECK(info == 0, sjoin("Failed the diagonalization of the updated GW 1-RDM with info:", itoa(info)))
391 :
392 : ! Sort in descending order
393 392 : do ib1dm=1,ndim
394 344 : occs_tmp2(ib1dm)=occs_tmp(ndim-(ib1dm-1))
395 2932 : do ib2dm=1,ndim
396 2932 : eigenvect(ib2dm,ib1dm)=tmp_mat(ib2dm,(ndim-(ib1dm-1)))
397 : end do
398 392 : if (abs(occs_tmp2(ib1dm))<tol10) then
399 106 : occs_tmp2(ib1dm)=zero
400 : end if
401 : end do
402 :
403 : ! Check orthonormality?
404 48 : if (check_Sijmat) then
405 0 : do ib1dm=1,ndim
406 0 : do ib2dm=1,ib1dm
407 : Sib1k_ib2k=czero
408 0 : do ib3dm=1,ndim
409 0 : Sib1k_ib2k=Sib1k_ib2k+conjg(eigenvect(ib3dm,ib1dm))*eigenvect(ib3dm,ib2dm)
410 : end do
411 0 : if (ib1dm==ib2dm) then
412 0 : if(abs(Sib1k_ib2k-cmplx(one,zero))>tol10) then
413 0 : write(msg,'(a45,i5,a1,i5,f10.5)') 'Large deviation from identity for bands ',ib1dm,' ',ib2dm,real(Sib1k_ib2k)
414 0 : call wrtout(std_out,msg)
415 : endif
416 : else
417 0 : if (abs(Sib1k_ib2k)>tol10) then
418 0 : write(msg,'(a45,i5,a1,i5,f10.5)') 'Large deviation from identity for bands ',ib1dm,' ',ib2dm,real(Sib1k_ib2k)
419 0 : call wrtout(std_out,msg)
420 : end if
421 : end if
422 : end do
423 : end do
424 : end if
425 :
426 : ! Print results
427 48 : if (info==0) then
428 48 : if (iinfo==0) then
429 27 : write(msg,'(a51,3f10.5)') 'Occs. after updating with Sx-Vxc corr. at k-point:',ebands%kptns(1:,ik_ibz)
430 : else
431 21 : write(msg,'(a51,3f10.5)') 'Occs. after updating with S_c correct. at k-point:',ebands%kptns(1:,ik_ibz)
432 : endif
433 48 : call wrtout(units, msg)
434 48 : ib1dm=ndim-(ndim/10)*10
435 48 : do ib2dm=1,(ndim/10)*10,10
436 0 : write(msg,'(f11.5,9f10.5)') occs_tmp2(ib2dm:ib2dm+9)
437 48 : call wrtout(units, msg)
438 : end do
439 48 : ib1dm=(ndim/10)*10+1
440 48 : write(msg,'(f11.5,*(f10.5))') occs_tmp2(ib1dm:)
441 48 : call wrtout(units, msg)
442 : else
443 0 : write(msg,'(a36,3f10.5)') 'Error computing occs. for k-point: ',ebands%kptns(1:,ik_ibz)
444 0 : call wrtout(units, msg)
445 : end if
446 :
447 : ! Store natural orbital eigenvectors matrix and occs. Also compute total number of electrons for this k-point
448 48 : toccs_k=zero
449 392 : do ib1dm=1,ndim
450 2932 : do ib2dm=1,ndim
451 2932 : nateigv(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=eigenvect(ib1dm,ib2dm)
452 : end do
453 344 : occs(ib1+(ib1dm-1),ik_ibz)=occs_tmp2(ib1dm) ! Overwrite the initial KS-DFT occs from ib1 to ib2
454 392 : toccs_k=toccs_k+occs_tmp2(ib1dm)
455 : end do
456 :
457 48 : write(msg,'(a22,i5,a3,i5,a21,f10.5)') ' Total occ. from band ',ib1,' to', ib2,' at current k-point: ',toccs_k
458 48 : call wrtout(units, msg)
459 48 : write(msg,'(a5)') ' '
460 48 : call wrtout(units, msg)
461 :
462 48 : ABI_FREE(rwork)
463 48 : ABI_FREE(work)
464 48 : ABI_FREE(tmp_mat)
465 48 : ABI_FREE(eigenvect)
466 48 : ABI_FREE(occs_tmp)
467 48 : ABI_FREE(occs_tmp2)
468 :
469 48 : end subroutine natoccs
470 : !!***
471 :
472 : !!****f* ABINIT/update_hdr_bst
473 : !! NAME
474 : !! update_hdr_bst
475 : !!
476 : !! FUNCTION
477 : !! Update the Hdr for the WFK and DEN files and the occ. numbers in the ebands file for a given k-point
478 : !!
479 : !! INPUTS
480 : !! Wfd<wfd_t>=Datatype gathering data on QP amplitudes.
481 : !! ngfft_in(18)=information on the fine FFT grid used for densities and potentials.
482 : !! b1gw=min band for given k in the interval where we update.
483 : !! b2gw=max band for given k in the interval where we update.
484 : !! occs= array containing the occ numbers for a given k-point occs_ks(nband,k-point).
485 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
486 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
487 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
488 : !!
489 : !! OUTPUT
490 : !! Updated Hdr and ebands information
491 : !!
492 : !! SOURCE
493 :
494 5 : subroutine update_hdr_bst(Wfd,occs,b1gw,b2gw,ebands,Hdr,ngfft_in)
495 :
496 : !Arguments ------------------------------------
497 : !scalars
498 : integer,intent(in) :: b1gw,b2gw
499 : integer,intent(in),dimension(3) :: ngfft_in
500 : type(ebands_t),target,intent(inout) :: ebands
501 : type(Hdr_type),intent(inout) :: Hdr
502 : class(wfd_t),intent(in) :: Wfd
503 : !arrays
504 : real(dp),intent(in) :: occs(:,:)
505 : !Local variables ------------------------------
506 : !scalars
507 : integer :: ib1dm,ib2dm,dim_bands,ikpoint
508 : !arrays
509 : !************************************************************************
510 :
511 : ! ebands occ (QP_ebands ones) are changed and never recoverd
512 35 : do ikpoint=1,ebands%nkpt
513 269 : ebands%occ(b1gw:b2gw,ikpoint,1) = occs(b1gw:b2gw,ikpoint) ! Spins summed, occ in [0:2]
514 : enddo
515 5 : ABI_COMMENT("QP_ebands: occupancies were updated with nat. orb. ones")
516 5 : if ((size(Hdr%occ(:))/ebands%nkpt) < (b2gw-b1gw+1)) then
517 : !Actually, we should never reach this point because the code should stop during Wfd initialization in m_sigma_driver
518 0 : ABI_ERROR("Impossible to use the existing read WFK to build a new one!")
519 : end if
520 :
521 : ! Update occ in Hdr before printing
522 5 : ib1dm=1
523 35 : do ikpoint=1,ebands%nkpt
524 30 : dim_bands=size(ebands%occ(:,ikpoint,1))
525 275 : do ib2dm=1,dim_bands
526 240 : Hdr%occ(ib1dm)=ebands%occ(ib2dm,ikpoint,1) ! Because Hdr%occ is a 1-D array
527 270 : ib1dm=ib1dm+1
528 : end do
529 : end do
530 :
531 35 : Hdr%npwarr(:)=Wfd%npwarr(:) ! Use the npw and ngfft = ones used in GW calc
532 20 : Hdr%ngfft(1:3)=ngfft_in(1:3)
533 5 : ABI_COMMENT("Hdr_sigma: occupancies, npw, and ngfft were updated")
534 :
535 5 : end subroutine update_hdr_bst
536 : !!***
537 :
538 : !!****f* ABINIT/print_tot_occ
539 : !! NAME
540 : !! print_tot_occ
541 : !!
542 : !! FUNCTION
543 : !! Compute and print the total (averaged) occ. from all k-points
544 : !!
545 : !! INPUTS
546 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
547 : !! eig(Sigp%nbnds,%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
548 : !! occ(Sigp%nbnds,%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
549 : !!
550 : !! OUTPUT
551 : !! Print the total (averaged) occ. = sum_k weight_k * Nelec_k
552 : !!
553 : !! SOURCE
554 :
555 5 : subroutine print_tot_occ(ebands)
556 :
557 : !Arguments ------------------------------------
558 : type(ebands_t),intent(in) :: ebands
559 :
560 : !Local variables-------------------------------
561 : !scalars
562 : character(len=500) :: msg
563 : integer :: ik,spin, units(2)
564 : real(dp) :: wtk,occ_bks,tot_occ
565 : ! *************************************************************************
566 :
567 15 : units = [std_out, ab_out]
568 :
569 5 : tot_occ=zero
570 :
571 10 : do spin=1,ebands%nsppol
572 40 : do ik=1,ebands%nkpt
573 30 : wtk = ebands%wtk(ik)
574 270 : occ_bks = sum(ebands%occ(:,ik,spin))
575 : !if (sigma%nsig_ab==1) then ! Only closed-shell restricted is programed
576 35 : tot_occ=tot_occ+occ_bks*wtk
577 : !end if
578 : end do
579 : end do
580 :
581 5 : write(msg,'(a1)') ' '
582 5 : call wrtout(units, msg)
583 5 : write(msg,'(a39,f10.5)') ' Total averaged occ. from all k-points: ',tot_occ
584 5 : call wrtout(units, msg)
585 5 : write(msg,'(a1)') ' '
586 5 : call wrtout(units, msg)
587 :
588 5 : end subroutine print_tot_occ
589 : !!***
590 :
591 : !!****f* ABINIT/get_chkprdm
592 : !! NAME
593 : !! get_chkprdm
594 : !!
595 : !! FUNCTION
596 : !! Read all checkpoint files built on previous runs
597 : !!
598 : !! INPUTS
599 : !! Wfd<wfd_t>=Wave function descriptor see file 69_wfd/m_wfd.F90
600 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
601 : !! Sigp<sigparams_t>=Parameters governing the self-energy calculation.
602 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
603 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
604 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
605 : !! occs = occ. numbers array occs(Wfd%mband,Wfd%nkibz)
606 : !! nateigv = natural orbital eigenvectors nateigv(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol))
607 : !! sigmak_todo = integer array initialized to 1 and its components are set to 0 if the kpoint
608 : !! is read from the checkpoint sigmak_todo(Wfd%nkibz)
609 : !! my_rank = rank of the mpi process.
610 : !!
611 : !! OUTPUT
612 : !! occ are updated if they are read from any checkpoint file
613 : !! nateigv are stored if they are read from any checkpoint file
614 : !! sigmak_todo components set to 1 if the kpoint is read from any checkpoint file
615 : !!
616 : !! SOURCE
617 :
618 1 : subroutine get_chkprdm(Wfd,Kmesh,Sigp,ebands,occs,nateigv,sigmak_todo,my_rank,gw1rdm_fname_in)
619 : !Arguments ------------------------------------
620 : !scalars
621 : integer,intent(in) :: my_rank
622 : class(wfd_t),intent(in) :: Wfd
623 : type(kmesh_t),intent(in) :: Kmesh
624 : type(sigparams_t),intent(in) :: Sigp
625 : type(ebands_t),intent(in) :: ebands
626 : character(len=fnlen),intent(in) :: gw1rdm_fname_in
627 : !arrays
628 : integer,intent(inout) :: sigmak_todo(:)
629 : real(dp),intent(inout) :: occs(:,:)
630 : complex(dp),intent(inout) :: nateigv(:,:,:,:)
631 : !Local variables-------------------------------
632 : !scalars
633 : integer,parameter :: master=0,iunit=666314
634 : integer :: ierr,ib1,ib2,ib3,ikcalc,istat,ik_ibz,ik_ibz_read,iread,iread_eigv
635 : real(dp) :: auxl_read
636 : character(len=fnlen) :: gw1rdm_fname
637 : character(len=500) :: msg
638 : !arrays
639 1 : real(dp),allocatable :: occ_tmp(:),eigvect_tmp(:)
640 : ! *************************************************************************
641 :
642 1 : if (my_rank==master) then
643 1 : iread_eigv=Wfd%mband
644 1 : iread_eigv=iread_eigv*(2*iread_eigv)
645 3 : ABI_MALLOC(occ_tmp,(Wfd%mband))
646 3 : ABI_MALLOC(eigvect_tmp,(iread_eigv))
647 :
648 7 : do ikcalc=1,Sigp%nkptgw
649 6 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
650 6 : if(ik_ibz<10) then
651 6 : write(gw1rdm_fname,"(a,i1)") trim(gw1rdm_fname_in),ik_ibz
652 0 : else if(ik_ibz<100 .and. ik_ibz>=10) then
653 0 : write(gw1rdm_fname,"(a,i2)") trim(gw1rdm_fname_in),ik_ibz
654 0 : else if(ik_ibz<1000 .and. ik_ibz>=100) then
655 0 : write(gw1rdm_fname,"(a,i3)") trim(gw1rdm_fname_in),ik_ibz
656 : else
657 0 : ABI_ERROR("The maximum k-point label for the checkpoint file to read is 999.")
658 : end if
659 6 : write(msg,'(a1)')' '
660 6 : call wrtout(std_out,msg)
661 6 : write(msg,'(a25,a)')' Reading checkpoint file ',trim(gw1rdm_fname)
662 6 : call wrtout(std_out,msg)
663 6 : write(msg,'(a1)')' '
664 6 : call wrtout(std_out,msg)
665 822 : occ_tmp(:)=zero;eigvect_tmp(:)=zero;
666 6 : open(unit=iunit,form='unformatted',file=gw1rdm_fname,iostat=istat,status='old')
667 6 : iread=0;ik_ibz_read=0;
668 6 : if (istat==0) then
669 : do
670 414 : if (iread<Wfd%mband) then
671 27 : iread=iread+1
672 27 : read(iunit,iostat=istat) auxl_read
673 27 : if (istat==0) then
674 24 : occ_tmp(iread)=auxl_read
675 : end if
676 387 : else if (iread<(iread_eigv+Wfd%mband)) then
677 384 : iread=iread+1
678 384 : read(iunit,iostat=istat) auxl_read
679 384 : if (istat==0) then
680 384 : eigvect_tmp(iread-Wfd%mband)=auxl_read
681 : end if
682 : else
683 3 : read(iunit,iostat=istat) ik_ibz_read
684 3 : if (istat==0 .and. ik_ibz_read/=0) then
685 3 : iread=0
686 3 : sigmak_todo(ik_ibz_read)=0
687 3 : ib3=1
688 27 : do ib1=1,Wfd%mband
689 24 : occs(ib1,ik_ibz_read)=occ_tmp(ib1)
690 219 : do ib2=1,Wfd%mband
691 192 : nateigv(ib1,ib2,ik_ibz_read,1)=cmplx(eigvect_tmp(ib3),eigvect_tmp(ib3+1))
692 216 : ib3=ib3+2
693 : end do
694 : end do
695 3 : ik_ibz_read=0
696 411 : occ_tmp=zero;eigvect_tmp=zero;
697 : end if
698 : end if
699 414 : if(istat/=0) then
700 : exit
701 : end if
702 : end do
703 : end if
704 7 : close(iunit)
705 : end do
706 1 : write(msg,'(a1)')' '
707 1 : call wrtout(std_out,msg)
708 1 : write(msg,'(a49)')' List of k-points read from all checkpoint files '
709 1 : call wrtout(std_out,msg)
710 7 : do ikcalc=1,Sigp%nkptgw
711 6 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Irred k-point for GW
712 7 : if (sigmak_todo(ik_ibz)==0) then
713 3 : write(msg,'(3f10.5)') ebands%kptns(1:,ik_ibz)
714 3 : call wrtout(std_out,msg)
715 : end if
716 : enddo
717 1 : write(msg,'(a1)')' '
718 1 : call wrtout(std_out,msg)
719 1 : ABI_FREE(occ_tmp)
720 1 : ABI_FREE(eigvect_tmp)
721 : end if
722 :
723 : ! Broadcast from master the information stored in occs and nateigv to all processes.
724 1 : call xmpi_barrier(Wfd%comm)
725 : ierr=0
726 1 : call xmpi_bcast(sigmak_todo(:),master,Wfd%comm,ierr)
727 1 : if(ierr/=0) then
728 0 : ABI_ERROR("Error distributing the sigmak_todo table.")
729 : endif
730 1 : call xmpi_bcast(occs(:,:),master,Wfd%comm,ierr)
731 1 : if(ierr/=0) then
732 0 : ABI_ERROR("Error distributing the occs read from checkpoint file(s).")
733 : endif
734 1 : call xmpi_bcast(nateigv(:,:,:,:),master,Wfd%comm,ierr)
735 1 : if(ierr/=0) then
736 0 : ABI_ERROR("Error distributing the natural orbital eigenvectors read from checkpoint file(s).")
737 : endif
738 :
739 1 : end subroutine get_chkprdm
740 : !!***
741 :
742 : !!****f* ABINIT/print_chkprdm
743 : !! NAME
744 : !! print_chkprdm
745 : !!
746 : !! FUNCTION
747 : !! Write the checkpoint file for a given k-point
748 : !!
749 : !! INPUTS
750 : !! Wfd<wfd_t>=Wave function descriptor see file 69_wfd/m_wfd.F90
751 : !! occs = occ. numbers array occs(Wfd%mband,Wfd%nkibz)
752 : !! nateigv = natural orbital eigenvectors nateigv(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol))
753 : !! my_rank = rank of the mpi process.
754 : !! gw1rdm_fname_out = name of the gw1rdm checkpoint out without k-point extension
755 : !!
756 : !! OUTPUT
757 : !!
758 : !! SOURCE
759 :
760 6 : subroutine print_chkprdm(Wfd,occs,nateigv,ik_ibz,my_rank,gw1rdm_fname_out)
761 : !Arguments ------------------------------------
762 : !scalars
763 : integer,intent(in) :: ik_ibz,my_rank
764 : class(wfd_t),intent(in) :: Wfd
765 : character(len=fnlen),intent(in) :: gw1rdm_fname_out
766 : !arrays
767 : real(dp),intent(in) :: occs(:,:)
768 : complex(dp),intent(in) :: nateigv(:,:,:,:)
769 : !Local variables-------------------------------
770 : !scalars
771 : integer,parameter :: master=0,iunit=666314
772 : integer :: iwrite,iwrite2
773 : character(len=fnlen) :: gw1rdm_fname
774 : character(len=500) :: msg
775 : !arrays
776 : ! *************************************************************************
777 :
778 6 : if (my_rank==master) then
779 6 : if(ik_ibz<10) then
780 6 : write(gw1rdm_fname,"(a,i1)") trim(gw1rdm_fname_out),ik_ibz
781 0 : else if(ik_ibz<100 .and. ik_ibz>=10) then
782 0 : write(gw1rdm_fname,"(a,i2)") trim(gw1rdm_fname_out),ik_ibz
783 0 : else if(ik_ibz<1000 .and. ik_ibz>=100) then
784 0 : write(gw1rdm_fname,"(a,i3)") trim(gw1rdm_fname_out),ik_ibz
785 : else
786 0 : ABI_ERROR("The maximum k-point label for the checkpoint file to write is 999.")
787 : end if
788 6 : write(msg,'(a1)')' '
789 6 : call wrtout(std_out,msg)
790 6 : write(msg,'(a25,a)')' Writing checkpoint file ',trim(gw1rdm_fname)
791 6 : call wrtout(std_out,msg)
792 6 : write(msg,'(a1)')' '
793 6 : call wrtout(std_out,msg)
794 6 : open(unit=iunit,form='unformatted',file=gw1rdm_fname)
795 54 : do iwrite=1,Wfd%mband
796 54 : write(iunit) occs(iwrite,ik_ibz)
797 : end do
798 54 : do iwrite=1,Wfd%mband
799 438 : do iwrite2=1,Wfd%mband
800 384 : write(iunit) real(nateigv(iwrite,iwrite2,ik_ibz,1))
801 432 : write(iunit) aimag(nateigv(iwrite,iwrite2,ik_ibz,1))
802 : end do
803 : end do
804 6 : write(iunit) ik_ibz
805 6 : close(iunit)
806 : end if
807 :
808 6 : call xmpi_barrier(Wfd%comm)
809 :
810 6 : end subroutine print_chkprdm
811 : !!***
812 :
813 : !!****f* ABINIT/change_matrix
814 : !! NAME
815 : !! change_matrix
816 : !!
817 : !! FUNCTION
818 : !! Transform integrals from KS -> NO and NO -> KS orbitals
819 : !!
820 : !! Transform <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j>,
821 : !! <KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j>,
822 : !! and <KS_i|T|KS_j> -> <NO_i|T|NO_j>
823 : !!
824 : !!
825 : !! INPUTS
826 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
827 : !! Sigp<sigparams_t>=Parameters governing the self-energy calculation.
828 : !! nateigv = natural orbital eigenvectors nateigv(Wfd%mband,Wfd%mband,Wfd%nkibz,Sigp%nsppol))
829 : !!
830 : !! OUTPUT
831 : !! Mels
832 : !! %kinetic=matrix elements of $t$.
833 : !! %vhartr =matrix elements of $v_H$.
834 : !! Sr=sigma_t (see the definition of this structured datatype)
835 : !!
836 : !! SOURCE
837 :
838 4 : subroutine change_matrix(Sigp,Sr,Mels,Kmesh,nateigv)
839 :
840 : !Arguments ------------------------------------
841 : !scalars
842 : type(kmesh_t),intent(in) :: Kmesh
843 : type(sigparams_t),intent(in) :: Sigp
844 : type(sigma_t),intent(inout) :: Sr
845 : type(melements_t),intent(inout) :: Mels
846 : !arrays
847 : complex(dp),intent(in) :: nateigv(:,:,:,:)
848 : !Local variables-------------------------------
849 : !scalars
850 : integer :: ikcalc,ik_ibz,ib1,ib2,ib1dm,ib2dm
851 : !arrays
852 4 : complex(dp),allocatable :: mat2rot(:,:),Umat(:,:)
853 : ! *************************************************************************
854 :
855 28 : do ikcalc=1,Sigp%nkptgw
856 24 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
857 48 : ib1=MINVAL(Sigp%minbnd(ikcalc,:)) ! min and max band indices for GW corrections (for this k-point)
858 48 : ib2=MAXVAL(Sigp%maxbnd(ikcalc,:))
859 96 : ABI_MALLOC(mat2rot,(ib2-ib1+1,ib2-ib1+1))
860 72 : ABI_MALLOC(Umat,(ib2-ib1+1,ib2-ib1+1))
861 : ! <NO_i|K[NO]|NO_j> -> <KS_i|K[NO]|KS_j>
862 216 : do ib1dm=1,ib2-ib1+1
863 1752 : do ib2dm=1,ib2-ib1+1
864 1536 : Umat(ib1dm,ib2dm)=nateigv(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)
865 1728 : mat2rot(ib1dm,ib2dm)=Sr%x_mat(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)
866 : end do
867 : end do
868 24 : call rotate_ks_no(ib1,ib2,mat2rot,Umat,0)
869 216 : do ib1dm=1,ib2-ib1+1
870 1752 : do ib2dm=1,ib2-ib1+1
871 1728 : Sr%x_mat(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=mat2rot(ib1dm,ib2dm)
872 : end do
873 : end do
874 : ! <KS_i|J[NO]|KS_j> -> <NO_i|J[NO]|NO_j>
875 216 : do ib1dm=1,ib2-ib1+1
876 1752 : do ib2dm=1,ib2-ib1+1
877 1728 : mat2rot(ib1dm,ib2dm)=Mels%vhartree(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)
878 : end do
879 : end do
880 24 : call rotate_ks_no(ib1,ib2,mat2rot,Umat,1)
881 216 : do ib1dm=1,ib2-ib1+1
882 1752 : do ib2dm=1,ib2-ib1+1
883 1728 : Mels%vhartree(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=mat2rot(ib1dm,ib2dm)
884 : end do
885 : end do
886 : ! <KS_i|T|KS_j> -> <NO_i|T|NO_j>
887 216 : do ib1dm=1,ib2-ib1+1
888 1752 : do ib2dm=1,ib2-ib1+1
889 1728 : mat2rot(ib1dm,ib2dm)=Mels%kinetic(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)
890 : end do
891 : end do
892 24 : call rotate_ks_no(ib1,ib2,mat2rot,Umat,1)
893 216 : do ib1dm=1,ib2-ib1+1
894 1752 : do ib2dm=1,ib2-ib1+1
895 1728 : Mels%kinetic(ib1+(ib1dm-1),ib1+(ib2dm-1),ik_ibz,1)=mat2rot(ib1dm,ib2dm)
896 : end do
897 : end do
898 24 : ABI_FREE(Umat)
899 28 : ABI_FREE(mat2rot)
900 : end do
901 4 : end subroutine change_matrix
902 : !!***
903 :
904 : !!****f* ABINIT/print_total_energy
905 : !! NAME
906 : !! print_total_energy
907 : !!
908 : !! FUNCTION
909 : !! Print total energy and energy components
910 : !!
911 : !! INPUTS
912 : !! all energy terms are self-explanatory
913 : !!
914 : !! OUTPUT
915 : !!
916 : !! SOURCE
917 :
918 4 : subroutine print_total_energy(ekin_energy,evext_energy,evextnl_energy,e_corepsp,eh_energy,ex_energy,&
919 : exc_mbb_energy,e_ewald,etot,etot2,den_int)
920 : !Arguments ------------------------------------
921 : !scalars
922 : real(dp),intent(in) :: ekin_energy,evext_energy,evextnl_energy,e_corepsp,eh_energy,ex_energy
923 : real(dp),intent(in) :: exc_mbb_energy,e_ewald,etot,etot2,den_int
924 :
925 : !Local variables-------------------------------
926 : character(len=500) :: msg
927 : integer :: units(2)
928 :
929 : !************************************************************************
930 :
931 12 : units = [std_out, ab_out]
932 :
933 4 : write(msg,'(a1)')' '
934 4 : call wrtout(units, msg)
935 : write(msg,'(a98)')'---------------------------------------------------------------&
936 4 : &----------------------------------'
937 4 : call wrtout(units, msg)
938 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Ekinetic = : ',ekin_energy,' Ha ,',ekin_energy*Ha_eV,' eV'
939 4 : call wrtout(units, msg)
940 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Evext_l = : ',evext_energy,' Ha ,',evext_energy*Ha_eV,' eV'
941 4 : call wrtout(units, msg)
942 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Evext_nl = : ',evextnl_energy,' Ha ,',evextnl_energy*Ha_eV,' eV'
943 4 : call wrtout(units, msg)
944 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Epsp_core = : ',e_corepsp,' Ha ,',e_corepsp*Ha_eV,' eV'
945 4 : call wrtout(units, msg)
946 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Ehartree = : ',eh_energy,' Ha ,',eh_energy*Ha_eV,' eV'
947 4 : call wrtout(units, msg)
948 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Ex[SD] = : ',ex_energy,' Ha ,',ex_energy*Ha_eV,' eV'
949 4 : call wrtout(units, msg)
950 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Exc[MBB] = : ',exc_mbb_energy,' Ha ,',exc_mbb_energy*Ha_eV,' eV'
951 4 : call wrtout(units, msg)
952 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Enn = : ',e_ewald,' Ha ,',e_ewald*Ha_eV,' eV'
953 4 : call wrtout(units, msg)
954 : write(msg,'(a98)')'-----------------------------------------------------------------&
955 4 : &--------------------------------'
956 4 : call wrtout(units, msg)
957 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Etot[SD] = : ',etot,' Ha ,',etot*Ha_eV,' eV'
958 4 : call wrtout(units, msg)
959 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Etot[MBB] = : ',etot2,' Ha ,',etot2*Ha_eV,' eV'
960 4 : call wrtout(units, msg)
961 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Vee[SD] = : ',(ex_energy+eh_energy),' Ha ,',(ex_energy+eh_energy)*Ha_eV,' eV'
962 4 : call wrtout(units, msg)
963 4 : write(msg,'(a,f10.5,a,f10.3,a)')' Vee[MBB] = : ',(exc_mbb_energy+eh_energy),' Ha ,',&
964 8 : &(exc_mbb_energy+eh_energy)*Ha_eV,' eV'
965 4 : call wrtout(units, msg)
966 4 : write(msg,'(a,f10.5)') ' Density = : ',den_int
967 4 : call wrtout(units, msg)
968 4 : write(msg,'(a)')' Vee[SD] (= Ehartree + Ex[SD]) energy obtained using GW 1-RDM:'
969 4 : call wrtout(units, msg)
970 4 : write(msg,'(a)')' Vee[MBB] (= Ehartree + Exc[MBB]) energy obtained using GW 1-RDM:'
971 4 : call wrtout(units, msg)
972 : write(msg,'(a98)')'-------------------------------------------------------------------&
973 4 : &------------------------------'
974 4 : call wrtout(units, msg)
975 :
976 4 : end subroutine print_total_energy
977 : !!***
978 :
979 : !!****f* ABINIT/print_band_energies
980 : !! NAME
981 : !! print_band_energies
982 : !!
983 : !! FUNCTION
984 : !! Print updated band energies
985 : !!
986 : !!
987 : !! INPUTS
988 : !! Kmesh <kmesh_t>=Structure describing the k-point sampling.
989 : !! Sigp<sigparams_t>=Parameters governing the self-energy calculation.
990 : !! Mels
991 : !! %kinetic=matrix elements of $t$.
992 : !! %vhartr =matrix elements of $v_H$.
993 : !! Sr=sigma_t (see the definition of this structured datatype)
994 : !! ebands=<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
995 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
996 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
997 : !!
998 : !! OUTPUT
999 : !!
1000 : !! SOURCE
1001 :
1002 4 : subroutine print_band_energies(b1gw,b2gw,Sr,Sigp,Mels,Kmesh,ebands,new_hartr,old_purex)
1003 : !Arguments ------------------------------------
1004 : !scalars
1005 : type(kmesh_t),intent(in) :: Kmesh
1006 : type(sigparams_t),intent(in) :: Sigp
1007 : type(sigma_t),intent(in) :: Sr
1008 : type(ebands_t),intent(in) :: ebands
1009 : type(melements_t),intent(in) :: Mels
1010 : integer,intent(in) :: b1gw,b2gw
1011 : !arrays
1012 : complex(dp),intent(in) :: old_purex(:,:),new_hartr(:,:)
1013 : !Local variables-------------------------------
1014 : !scalars
1015 : integer :: ib,ikcalc,ik_ibz, units(2)
1016 : real(dp) :: eik_new
1017 : complex(dp) :: delta_band_ibik
1018 : character(len=500) :: msg
1019 : !************************************************************************
1020 :
1021 12 : units = [std_out, ab_out]
1022 :
1023 4 : write(msg,'(a1)') ' '
1024 4 : call wrtout(units, msg)
1025 4 : write(msg,'(a42)') ' Computing band corrections Delta eik (eV)'
1026 4 : call wrtout(units, msg)
1027 4 : write(msg,'(a42)') ' -----------------------------------------'
1028 4 : call wrtout(units, msg)
1029 4 : write(msg,'(a1)') ' '
1030 4 : call wrtout(units, msg)
1031 4 : write(msg,'(a1)') ' '
1032 4 : call wrtout(units, msg)
1033 : write(msg,'(a110)') ' Band corrections Delta eik = <KS_i|K[NO]-a*K[KS]+vH[NO]&
1034 4 : &-vH[KS]-Vxc[KS]|KS_i> and eik^new = eik^GS + Delta eik'
1035 4 : call wrtout(units, msg)
1036 4 : write(msg,'(a1)') ' '
1037 4 : call wrtout(units, msg)
1038 28 : do ikcalc=1,Sigp%nkptgw
1039 24 : ik_ibz=Kmesh%tab(Sigp%kptgw2bz(ikcalc)) ! Index of the irreducible k-point for GW
1040 : write(msg,'(a127)')'---------------------------------------------------------&
1041 24 : &--------------------------------------------------------------------'
1042 24 : call wrtout(units, msg)
1043 : write(msg,'(a)')' k-point band eik^GS eik^new Delta eik &
1044 24 : & K[NO] a*K[KS] Vxc[KS] vH[NO] vH[KS]'
1045 24 : call wrtout(units, msg)
1046 220 : do ib=b1gw,b2gw
1047 : delta_band_ibik=(new_hartr(ib,ikcalc)-Mels%vhartree(ib,ib,ik_ibz,1))&
1048 192 : &+Sr%x_mat(ib,ib,ik_ibz,1)-Mels%vxcval(ib,ib,ik_ibz,1)-old_purex(ib,ikcalc)
1049 192 : eik_new=real(ebands%eig(ib,ik_ibz,1))+real(delta_band_ibik)
1050 : write(msg,'(i5,4x,i5,8(4x,f10.3))') &
1051 192 : & ik_ibz,ib,real(ebands%eig(ib,ik_ibz,1))*Ha_eV,eik_new*Ha_eV,real(delta_band_ibik)*Ha_eV,&
1052 192 : & real(Sr%x_mat(ib,ib,ik_ibz,1))*Ha_eV,real(old_purex(ib,ikcalc))*Ha_eV,&
1053 192 : & real(Mels%vxcval(ib,ib,ik_ibz,1))*Ha_eV,&
1054 384 : & real(new_hartr(ib,ikcalc))*Ha_eV,real(Mels%vhartree(ib,ib,ik_ibz,1))*Ha_eV
1055 216 : call wrtout(units, msg)
1056 : enddo
1057 : enddo
1058 : write(msg,'(a127)')'---------------------------------------------------------&
1059 4 : &--------------------------------------------------------------------'
1060 4 : call wrtout(units, msg)
1061 :
1062 4 : end subroutine print_band_energies
1063 : !!***
1064 :
1065 : !!****f* ABINIT/rotate_ks_no
1066 : !! NAME
1067 : !! rotate_ks_no
1068 : !!
1069 : !! FUNCTION
1070 : !! Rotate a matrix from KS to NO basis and vicerversa.
1071 : !!
1072 : !! INPUTS
1073 : !! ib1=min band for given k
1074 : !! ib2=max band for given k.
1075 : !! Umat=array containing the eigenvectors in Columns (a unitary matrix)
1076 : !! Mat=initially an array containing the matrix elements in KS or NO basis
1077 : !! option=0 rotate from NO -> KS | 1 rotate from KS -> NO
1078 : !!
1079 : !! OUTPUT
1080 : !! Rotate a matrix from KS to NO basis and vicerversa and save the new matrix on Mat.
1081 : !! Mat=at the end an array containing the matrix elements in NO or KS basis
1082 : !!
1083 : !! SOURCE
1084 :
1085 72 : subroutine rotate_ks_no(ib1,ib2,Mat,Umat,option)
1086 : !Arguments ------------------------------------
1087 : !scalars
1088 : integer,intent(in) :: ib1,ib2,option
1089 : !arrays
1090 : complex(dp),intent(in) :: Umat(:,:)
1091 : complex(dp),intent(inout) :: Mat(:,:)
1092 : !Local variables ------------------------------
1093 : !scalars
1094 : integer:: ndim
1095 : !arrays
1096 : !************************************************************************
1097 :
1098 72 : ndim=ib2-ib1+1
1099 72 : if (option==0) then
1100 24 : call no2ks(ndim,Mat,Umat)
1101 : else
1102 48 : call ks2no(ndim,Mat,Umat)
1103 : end if
1104 :
1105 72 : end subroutine rotate_ks_no
1106 : !!***
1107 :
1108 : !!****f* ABINIT/ks2no
1109 : !! NAME
1110 : !! ks2no
1111 : !!
1112 : !! FUNCTION
1113 : !! Transform the matrix mat from KS to NO basis
1114 : !!
1115 : !! INPUTS
1116 : !! dim=dimension of the matrices
1117 : !! mat=array in the KS basis
1118 : !! rot=unitary matrix containg the eigenvectors in NO basis
1119 : !!
1120 : !! OUTPUT
1121 : !! mat=array in the NO basis
1122 : !!
1123 : !! SOURCE
1124 :
1125 48 : subroutine ks2no(ndim,mat,rot)
1126 : !Arguments ------------------------------------
1127 : !scalars
1128 : integer,intent(in) :: ndim
1129 : !arrays
1130 : complex(dp),dimension(:,:),intent(in) :: rot
1131 : complex(dp),dimension(:,:),intent(inout) :: mat
1132 : !Local variables ------------------------------
1133 : !scalars
1134 : !arrays
1135 48 : complex(dp),allocatable :: res(:,:)
1136 : !************************************************************************
1137 :
1138 192 : ABI_MALLOC(res,(ndim,ndim))
1139 3504 : res=czero
1140 :
1141 : ! <NO|Op|NO> = (U^t)* <KS|Op|KS> U
1142 31632 : res=matmul(conjg(transpose(rot)),mat)
1143 56208 : mat=matmul(res,rot)
1144 :
1145 48 : ABI_FREE(res)
1146 :
1147 48 : end subroutine ks2no
1148 : !!***
1149 :
1150 : !!****f* ABINIT/no2ks
1151 : !! NAME
1152 : !! no2ks
1153 : !!
1154 : !! FUNCTION
1155 : !! Transform the matrix mat from NO to KS basis
1156 : !!
1157 : !! INPUTS
1158 : !! dim=dimension of the matrices
1159 : !! mat=array in the KS basis
1160 : !! rot=unitary matrix containg the eigenvectors in NO basis
1161 : !!
1162 : !! OUTPUT
1163 : !! mat=array in the KS basis
1164 : !!
1165 : !! SOURCE
1166 :
1167 24 : subroutine no2ks(ndim,mat,rot)
1168 : !Arguments ------------------------------------
1169 : !scalars
1170 : integer,intent(in) :: ndim
1171 : !arrays
1172 : complex(dp),dimension(:,:),intent(in) :: rot
1173 : complex(dp),dimension(:,:),intent(inout) :: mat
1174 : !Local variables ------------------------------
1175 24 : complex(dp),allocatable :: res(:,:)
1176 : !************************************************************************
1177 :
1178 96 : ABI_MALLOC(res,(ndim,ndim))
1179 1752 : res=czero
1180 :
1181 : ! <KS|Op|KS> = U <NO|Op|NO> (U^t)*
1182 15816 : res=matmul(rot,mat)
1183 28104 : mat=matmul(res,conjg(transpose(rot)))
1184 :
1185 24 : ABI_FREE(res)
1186 :
1187 24 : end subroutine no2ks
1188 : !!***
1189 :
1190 : !!****f* ABINIT/printrdm_k
1191 : !! NAME
1192 : !! printrdm_k
1193 : !!
1194 : !! FUNCTION
1195 : !! Print the DM1 matrix
1196 : !!
1197 : !! INPUTS
1198 : !! ib1=min band.
1199 : !! ib2=max band.
1200 : !! rdm_k=array containing the 1-RDM matrix
1201 : !!
1202 : !! OUTPUT
1203 : !! Print the 1-RDM matrix
1204 : !! SOURCE
1205 :
1206 : subroutine printrdm_k(ib1,ib2,rdm_k) ! Only used for debug on this file, do not use it with large arrays!
1207 : !Arguments ------------------------------------
1208 : !scalars
1209 : integer,intent(in) :: ib1,ib2
1210 : !arrays
1211 : complex(dp),intent(in) :: rdm_k(:,:)
1212 : !Local variables ------------------------------
1213 : !scalars
1214 : integer::ib1dm
1215 : character(len=500) :: msg
1216 : !arrays
1217 : !************************************************************************
1218 :
1219 : do ib1dm=ib1,ib2
1220 : write(msg,'(*(f12.5))') real(rdm_k(ib1dm,ib1:ib2))
1221 : call wrtout(std_out, msg)
1222 : end do
1223 :
1224 : end subroutine printrdm_k
1225 : !!***
1226 :
1227 : end module m_gwrdm
1228 : !!***
|
