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