Line data Source code
1 : !!****m* ABINIT/m_datafordmft
2 : !! NAME
3 : !! m_datafordmft
4 : !!
5 : !! FUNCTION
6 : !! This module produces inputs for the DMFT calculation
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! INPUTS
15 : !!
16 : !! OUTPUT
17 : !!
18 : !! SOURCE
19 :
20 : #if defined HAVE_CONFIG_H
21 : #include "config.h"
22 : #endif
23 :
24 : #include "abi_common.h"
25 :
26 : MODULE m_datafordmft
27 :
28 : use defs_abitypes
29 : use defs_basis
30 : use m_abi_linalg, only : abi_xgemm
31 : use m_crystal, only : crystal_t
32 : use m_dtset
33 : use m_errors
34 : use m_fstrings, only : int2char4
35 : use m_io_tools, only : open_file
36 : use m_matlu, only : add_matlu,checkdiag_matlu,destroy_matlu,diff_matlu,copy_matlu, &
37 : & init_matlu,matlu_type,print_matlu,sym_matlu,xmpi_matlu,magnfield_matlu
38 : use m_matrix, only : invsqrt_matrix
39 : use m_mpinfo, only : proc_distrb_cycle
40 : use m_oper, only : copy_oper,destroy_oper,diff_oper,downfold_oper,identity_oper,init_oper,oper_type,prod_oper
41 : use m_paw_dmft, only : paw_dmft_type
42 : use m_paw_ij, only : paw_ij_type
43 : use m_pawcprj, only : pawcprj_alloc,pawcprj_free,pawcprj_get,pawcprj_type
44 : use m_pawtab, only : pawtab_type
45 : use m_xmpi
46 :
47 : implicit none
48 :
49 : private
50 :
51 : public :: datafordmft
52 : public :: chipsi_print
53 : public :: compute_levels
54 : public :: chipsi_renormalization
55 : public :: hybridization_asymptotic_coefficient
56 : public :: compute_wannier
57 : public :: print_wannier
58 : !!***
59 :
60 : contains
61 :
62 : !!****f* ABINIT/datafordmft
63 : !! NAME
64 : !! datafordmft
65 : !!
66 : !! FUNCTION
67 : !! Compute chipsi (and print some data for check)
68 : !!
69 : !! INPUTS
70 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions.
71 : !! cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
72 : !! and each |p_lmn> non-local projector
73 : !! cryst_struc <type(crystal_t)>=crystal structure data
74 : !! -gprimd(3,3)=dimensional reciprocal space primitive translations
75 : !! -indsym(4,nsym,natom)=indirect indexing array for atom labels
76 : !! -symrec(3,3,nsym)=symmetry operations in reciprocal space
77 : !! - nsym= number of symetry operations
78 : !! dft_occup <type(oper_type)> = DFT occupations of the correlated orbitals
79 : !! dimcprj(natom) = dimension for cprj
80 : !! dtset <type(dataset_type)>=all input variables for this dataset
81 : !! eigen(mband*nkpt*nsppol)=array for holding DFT eigenvalues (hartree)
82 : !! mband_cprj=number of bands on each process of the band communicator
83 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
84 : !! mpi_enreg=information about MPI parallelization
85 : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
86 : !! occ(mband*nkpt*nsppol) = occupancies of KS states.
87 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
88 : !! paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
89 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
90 : !! usecprj=1 if cprj datastructure is stored in memory
91 : !! nbandkss=number of bands in the KSS file
92 : !!
93 : !! OUTPUT
94 : !! paw_dmft%chipsi((2*maxlpawu+1)*nspinor,mbandc,nkpt,nsppol,natom): projections <Chi|Psi>
95 : !! paw_dmft%eigen(paw_dmft%mbandc,paw_dmft%nkpt,paw_dmft%nsppol)
96 : !!
97 : !! SIDE EFFECTS
98 : !! (only writing, printing)
99 : !!
100 : !! NOTES
101 : !!
102 : !! SOURCE
103 :
104 97 : subroutine datafordmft(cg,cprj,cryst_struc,dft_occup,dimcprj,dtset,eigen,mband_cprj,mcg,&
105 97 : & mpi_enreg,my_nspinor,occ,paw_dmft,paw_ij,pawtab,usecprj,nbandkss)
106 :
107 : !Arguments ------------------------------------
108 : integer, intent(in) :: mband_cprj,mcg,my_nspinor,usecprj
109 : integer, optional, intent(in) :: nbandkss
110 : type(MPI_type), intent(in) :: mpi_enreg
111 : type(dataset_type), intent(in) :: dtset
112 : type(oper_type), intent(inout) :: dft_occup !vz_i
113 : type(crystal_t), intent(in) :: cryst_struc
114 : type(paw_dmft_type), intent(inout) :: paw_dmft
115 : integer, intent(in) :: dimcprj(paw_dmft%natom)
116 : real(dp), intent(in) :: occ(paw_dmft%mband*paw_dmft%nkpt*paw_dmft%nsppol)
117 : real(dp), target, intent(in) :: eigen(paw_dmft%mband*paw_dmft%nkpt*paw_dmft%nsppol)
118 : real(dp), intent(in) :: cg(2,mcg)
119 : type(paw_ij_type), intent(in) :: paw_ij(paw_dmft%natom)
120 : ! type(pawcprj_type) :: cprj(cryst_struc%natom,my_nspinor*mband*mkmem*nsppol)
121 : type(pawcprj_type), intent(in) :: cprj(paw_dmft%natom,my_nspinor*mband_cprj*dtset%mkmem*paw_dmft%nsppol*usecprj)
122 : type(pawtab_type), intent(in) :: pawtab(paw_dmft%ntypat)
123 : !Local variables-------------------------------
124 : integer :: band_index,comm_band,comm_kpt,iatom,ib,iband,ibandc,ibuf_chipsi,ibuf_psi,icg,icgb
125 : integer :: icprj,idijeff,ierr,ik,ikpt,ilmn,im,iorder_cprj,iproj,ir,irank,ispinor,ispinor1
126 : integer :: isppol,itypat,lmn_size,lpawu,lpawu1,maxlpawu,maxmeshsize,maxnproju,mband,mbandc
127 : integer :: me_band,me_kpt,mkmem,natom,nband_k,nband_k_cprj,nbandf,nbandi,ndim,nkpt
128 : integer :: nproc_band,nproc_spkpt,nproju,npw,nspinor,nsploop,nsppol,nsppol_mem,opt_renorm
129 : integer :: option,paral_kgb,pawprtvol,siz_buf,siz_buf_psi,siz_paw,siz_proj,siz_wan,unt
130 : logical :: prt_wan,t2g,use_full_chipsi,verif_band,x2my2d
131 : real(dp) :: rint,bfield
132 : character(len=500) :: message
133 97 : type(oper_type) :: loc_norm_check
134 97 : integer, allocatable :: displs(:),recvcounts(:)
135 97 : complex(dp), allocatable :: buf_chipsi(:),buf_chipsi_tot(:),chipsi_tmp(:),cwprj(:,:)
136 97 : type(pawcprj_type), allocatable :: cwaveprj(:,:)
137 97 : type(matlu_type), allocatable :: matlu_temp(:)
138 : integer, parameter :: spinor_idxs(2,4) = RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
139 : !************************************************************************
140 :
141 : !DBG_ENTER("COLL")
142 :
143 97 : mband = paw_dmft%mband
144 97 : mbandc = paw_dmft%mbandc
145 97 : mkmem = paw_dmft%mkmem
146 97 : natom = paw_dmft%natom
147 97 : nkpt = paw_dmft%nkpt
148 97 : nspinor = paw_dmft%nspinor
149 97 : nsppol = paw_dmft%nsppol
150 232 : nsppol_mem = sum(mpi_enreg%my_isppoltab(1:nsppol))
151 97 : pawprtvol = dtset%pawprtvol
152 97 : prt_wan = (paw_dmft%dmft_prtwan == 1)
153 :
154 97 : if (abs(pawprtvol) >= 3) then
155 73 : write(message,*) " number of k-points used is nkpt=nkpt ",nkpt
156 73 : call wrtout(std_out,message,'COLL')
157 73 : write(message,*) " warning: parallelized version ",nkpt
158 73 : call wrtout(std_out,message,'COLL')
159 73 : write(message,*) " weights k-points used is wtk=wtk"
160 73 : call wrtout(std_out,message,'COLL')
161 : end if ! abs(pawprtvol)>=3
162 :
163 97 : if (usecprj == 0) then
164 0 : write(message,*) " usecprj=0 : BUG in datafordmft",usecprj
165 0 : ABI_BUG(message)
166 : end if
167 :
168 97 : if (my_nspinor /= nspinor) then
169 0 : write(message,*) " my_nspinor=/dtset%nspinor, datafordmft not working in this case",my_nspinor,nspinor
170 0 : ABI_ERROR(message)
171 : end if
172 :
173 :
174 : !do ib=1,my_nspinor*mband_cprj*mkmem*nsppol*usecprj
175 : !write(std_out,'(a,i6,3e16.7)') "cprj",ib,cprj(1,ib)%cp(1,19),cprj(1,ib)%cp(2,19),cprj(1,ib)%cp(1,19)**2+cprj(1,ib)%cp(2,19)**2
176 : !enddo
177 :
178 : !----------------------------------- MPI-------------------------------------
179 :
180 : ! Init parallelism
181 97 : paral_kgb = mpi_enreg%paral_kgb
182 97 : comm_kpt = merge(mpi_enreg%comm_kpt,mpi_enreg%comm_cell,paral_kgb==1)
183 97 : comm_band = mpi_enreg%comm_band
184 97 : me_kpt = mpi_enreg%me_kpt
185 97 : me_band = mpi_enreg%me_band
186 97 : nproc_band = mpi_enreg%nproc_band
187 97 : nproc_spkpt = mpi_enreg%nproc_spkpt
188 :
189 97 : if (nproc_band /= (mband/mband_cprj)) then
190 0 : message = "Inconsistency in datafordmft: nproc_band should be equal to mband/mband_cprj"
191 0 : ABI_BUG(message)
192 : end if
193 :
194 97 : iorder_cprj = 0
195 97 : ABI_CHECK(dtset%mkmem/=0,"mkmem=0 not supported anymore!")
196 : !todo_ab: extract cprj from file unpaw in the following..
197 : !call abi_abort('COLL')
198 :
199 : !----------------------------------- MPI-------------------------------------
200 :
201 97 : nbandi = paw_dmft%dmftbandi
202 97 : nbandf = paw_dmft%dmftbandf
203 97 : t2g = (paw_dmft%dmft_t2g == 1)
204 97 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
205 : !use_full_chipsi = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
206 97 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
207 : use_full_chipsi = .true.
208 : else
209 97 : use_full_chipsi = (paw_dmft%dmft_full_chipsi /= 0)
210 : endif
211 :
212 97 : if (use_full_chipsi .and. mpi_enreg%nproc_fft > 1) then
213 0 : message = "datafordmft not working when nproc_fft > 1 and use_full_chipsi=1"
214 0 : ABI_ERROR(message)
215 : end if
216 :
217 : !if(mpi_enreg%me==0) write(7886,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc
218 : !if(mpi_enreg%me==1) write(7887,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc
219 : !if(mpi_enreg%me==2) write(7888,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc
220 97 : write(message,'(2a)') ch10,' == Prepare data for DFT+DMFT calculation'
221 97 : call wrtout(std_out,message,'COLL')
222 97 : if (abs(pawprtvol) >= 3) then
223 73 : write(message,'(2a)') ch10,'---------------------------------------------------------------'
224 : ! call wrtout(ab_out,message,'COLL')
225 73 : call wrtout(std_out,message,'COLL')
226 73 : write(message,'(12a)') ch10,' Print useful data (as a check)',ch10,&
227 73 : & ' - Overlap of KS wfc with atomic orbital inside sphere',ch10,&
228 73 : & ' - Eigenvalues',ch10,&
229 73 : & ' - Weights of k-points',ch10,&
230 73 : & ' - Number of spins ',ch10,&
231 146 : & ' - Number of states'
232 : ! call wrtout(ab_out,message,'COLL')
233 73 : call wrtout(std_out,message,'COLL')
234 73 : write(message,'(2a)') ch10,'---------------------------------------------------------------'
235 73 : call wrtout(std_out,message,'COLL')
236 : end if ! abs(pawprtvol)>=3
237 :
238 97 : if (dtset%nstep == 0 .and. dtset%nbandkss == 0) then
239 0 : message = 'nstep should be greater than 1'
240 0 : ABI_BUG(message)
241 : end if
242 :
243 : !********************* Max Values for U terms.
244 97 : maxlpawu = paw_dmft%maxlpawu
245 97 : maxmeshsize = paw_dmft%maxmeshsize
246 97 : maxnproju = paw_dmft%maxnproju
247 :
248 : !***************** in forlb.eig
249 97 : if (paw_dmft%myproc == 0 .and. abs(pawprtvol) >= 3) then
250 25 : if (open_file('forlb.eig',message,newunit=unt,form='formatted',status='unknown') /= 0) ABI_ERROR(message)
251 25 : rewind(unt)
252 25 : write(unt,*) "Number of bands, spins, and k-point; and spin-orbit flag"
253 25 : write(unt,*) mband,nsppol,nkpt,my_nspinor,nbandi,nbandf
254 25 : write(unt,*) " For each k-point, eigenvalues for each band"
255 258 : write(unt,*) (dtset%wtk(ikpt),ikpt=1,nkpt)
256 25 : band_index = 0
257 55 : do isppol=1,nsppol
258 30 : write(unt,*) " For spin"
259 30 : write(unt,*) isppol
260 378 : do ikpt=1,nkpt
261 323 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
262 323 : ibandc = 0
263 323 : write(unt,*) " For k-point"
264 323 : write(unt,*) ikpt
265 7427 : do iband=1,mband
266 7427 : if (paw_dmft%band_in(iband)) then
267 3003 : ibandc = ibandc + 1
268 3003 : write(unt,'(2i6,4x,f20.15)') ibandc,ikpt,eigen(iband+band_index)*two
269 : end if
270 : end do ! iband
271 353 : band_index = band_index + nband_k
272 : end do ! ikpt
273 : end do ! isppol
274 25 : close(unt)
275 : end if ! proc=me
276 :
277 : !== Setup Zeeman Contributios -mu_b.ge.Sz.Bz to Kohn-Sham energies for nspinor==1
278 97 : if(paw_dmft%dmft_magnfield .eq. 1) then
279 9 : if(nspinor .eq. 2) then
280 0 : write(message,'(a,a,2a)') ch10,&
281 0 : & 'Applying Zeeman contributions to Kohn-Sham energies is only implemented for nspinor = 1.'
282 0 : ABI_ERROR(message)
283 : else
284 9 : bfield=paw_dmft%dmft_magnfield_b
285 9 : write(message,'(2a)') ch10,'Adding Zeeman contribution to DFT eigenvalues'
286 9 : call wrtout(std_out,message,'COLL')
287 : endif
288 : else
289 : bfield=0.0
290 : endif
291 :
292 : !== put eigen into eigen_dft
293 97 : paw_dmft%eigen => eigen(:)
294 97 : band_index = 0
295 232 : do isppol=1,nsppol
296 1354 : do ikpt=1,nkpt
297 1122 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
298 1122 : ibandc = 0
299 27042 : do iband=1,mband
300 27042 : if (paw_dmft%band_in(iband)) then
301 10248 : ibandc = ibandc + 1
302 : ! paw_dmft%eigen_dft(ibandc,ikpt,isppol) = eigen(iband+band_index) ! in Ha
303 : ! paw_dmft%eigen_dft(isppol,ikpt,ibandc)=fermie
304 10248 : if (isppol .eq. 1) then
305 7554 : paw_dmft%eigen_dft(ibandc,ikpt,isppol) = eigen(iband+band_index)-half*bfield ! in Ha
306 2694 : elseif(isppol .eq. 2) then
307 2694 : paw_dmft%eigen_dft(ibandc,ikpt,isppol) = eigen(iband+band_index)+half*bfield ! in Ha
308 : endif
309 : end if
310 : end do ! iband
311 1257 : band_index = band_index + nband_k
312 : end do ! ikpt
313 : end do ! isppol
314 :
315 97 : if (abs(pawprtvol) >= 3) then
316 73 : write(message,'(2a)') ch10,' datafordmft : eigenvalues written on file'
317 73 : call wrtout(std_out,message,'COLL')
318 : end if
319 : !==========================================================================
320 : !***************** Compute <Chi|Psi>= <Chi|Psi_tilde> + \sum_{proja} <P_a|Psi><Chi|phi_a-phi_tilde_a>
321 : !==========================================================================
322 : !write(std_out,*) "size(cprj,dim=1)",size(cprj,dim=1),size(cprj,dim=2),dtset%mband,dtset%mkmem,dtset%nkpt
323 :
324 : !Allocate temporary cwaveprj storage
325 812 : ABI_MALLOC(cwaveprj,(natom,my_nspinor))
326 : !write(std_out,*) "before alloc cprj"
327 : !write(std_out,*) size(cwaveprj,dim=1),size(cwaveprj,dim=2),size(dimcprj,dim=1)
328 :
329 97 : call pawcprj_alloc(cwaveprj(:,:),0,dimcprj(:))
330 : !write(std_out,*) "after alloc cprj"
331 :
332 388 : ABI_MALLOC(cwprj,(maxnproju,2*maxlpawu+1))
333 291 : ABI_MALLOC(chipsi_tmp,(maxmeshsize))
334 :
335 97 : siz_buf = 0
336 97 : siz_buf_psi = 0
337 422 : do iatom=1,natom
338 325 : lpawu = paw_dmft%lpawu(iatom)
339 325 : if (lpawu == -1) cycle
340 109 : itypat = paw_dmft%typat(iatom)
341 109 : ndim = 2*lpawu + 1
342 109 : siz_buf = siz_buf + ndim
343 422 : siz_buf_psi = siz_buf_psi + ndim*paw_dmft%radgrid(itypat)%mesh_size
344 : end do ! iatom
345 97 : siz_buf = siz_buf * mbandc * nspinor * mkmem * nsppol_mem
346 97 : siz_buf_psi = siz_buf_psi * mbandc * nspinor * mkmem * nsppol_mem
347 :
348 291 : ABI_MALLOC(recvcounts,(nproc_spkpt))
349 194 : ABI_MALLOC(displs,(nproc_spkpt))
350 97 : call xmpi_allgather(siz_buf,recvcounts(:),comm_kpt,ierr)
351 :
352 97 : displs(1) = 0
353 193 : do irank=2,nproc_spkpt
354 193 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
355 : end do ! irank
356 :
357 291 : ABI_MALLOC(buf_chipsi,(siz_buf))
358 291 : ABI_MALLOC(buf_chipsi_tot,(recvcounts(nproc_spkpt)+displs(nproc_spkpt)))
359 97 : if (prt_wan) then
360 0 : ABI_MALLOC(paw_dmft%buf_psi,(siz_buf_psi))
361 : end if
362 :
363 97 : icprj = 0
364 97 : icg = 0
365 97 : ibuf_psi = 0
366 97 : ibuf_chipsi = 0
367 :
368 31568 : buf_chipsi(:) = czero
369 :
370 232 : do isppol=1,nsppol
371 :
372 135 : if (mpi_enreg%my_isppoltab(isppol) == 0) cycle
373 115 : ik = 0
374 :
375 1070 : do ikpt=1,nkpt
376 :
377 858 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
378 858 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle
379 :
380 633 : nband_k_cprj = nband_k / nproc_band
381 633 : npw = paw_dmft%npwarr(ikpt)
382 633 : icgb = icg
383 633 : ik = ik + 1
384 633 : ibandc = 0
385 : ! LOOP OVER BANDS
386 633 : ib = 0
387 :
388 16449 : do iband=1,nband_k
389 :
390 : ! Parallelization: treat only some bands
391 15816 : verif_band = .true.
392 15816 : if (paral_kgb == 1) then
393 6400 : if (mod((iband-1)/mpi_enreg%bandpp,nproc_band) /= me_band) verif_band = .false.
394 : else
395 9416 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol) /= me_kpt) verif_band = .false.
396 : end if
397 :
398 11016 : if (verif_band) ib = ib + 1
399 :
400 15816 : if (paw_dmft%band_in(iband)) then
401 6075 : ibandc = ibandc + 1
402 : else
403 9741 : icgb = icgb + npw*nspinor
404 9741 : cycle
405 : end if
406 :
407 6075 : if (verif_band) then
408 : call pawcprj_get(cryst_struc%atindx1(:),cwaveprj(:,:),cprj(:,:),natom,ib,icprj,ikpt,&
409 : & iorder_cprj,isppol,mband_cprj,dtset%mkmem,natom,1,nband_k_cprj,&
410 : & my_nspinor,nsppol,paw_dmft%unpaw,mpicomm=mpi_enreg%comm_kpt,&
411 4275 : & proc_distrb=mpi_enreg%proc_distrb(:,:,:))
412 : end if
413 :
414 12847 : do ispinor=1,my_nspinor
415 25814 : do iatom=1,natom
416 19675 : lpawu = paw_dmft%lpawu(iatom)
417 19675 : if (lpawu == -1) cycle
418 6851 : lpawu1 = merge(2,lpawu,t2g.or.x2my2d)
419 6851 : itypat = paw_dmft%typat(iatom)
420 6851 : ndim = 2*lpawu + 1
421 6851 : nproju = pawtab(itypat)%nproju
422 6851 : siz_proj = paw_dmft%siz_proj(itypat)
423 6851 : siz_wan = paw_dmft%radgrid(itypat)%mesh_size
424 6851 : siz_paw = min(siz_wan,paw_dmft%int_meshsz(itypat))
425 6851 : rint = paw_dmft%radgrid(itypat)%rad(siz_proj)
426 :
427 6851 : if (verif_band) then
428 5051 : lmn_size = pawtab(itypat)%lmn_size
429 112321 : do ilmn=1,lmn_size
430 : ! ------------ Select l=lpawu.
431 107270 : if (pawtab(itypat)%indlmn(1,ilmn) /= lpawu1) cycle
432 :
433 48502 : im = pawtab(itypat)%indlmn(2,ilmn) + lpawu1 + 1
434 48502 : if (x2my2d) then
435 0 : if (im /= 5) cycle
436 : im = 1
437 48502 : else if (t2g) then
438 10200 : if (im == 3 .or. im == 5) cycle
439 6120 : if (im == 4) im = 3
440 : end if
441 44422 : iproj = pawtab(itypat)%indlmn(3,ilmn)
442 : cwprj(iproj,im) = cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn), &
443 112321 : & cwaveprj(iatom,ispinor)%cp(2,ilmn),kind=dp)
444 : end do ! ilmn
445 : end if ! verif
446 :
447 44461 : do im=1,ndim
448 :
449 31471 : ibuf_chipsi = ibuf_chipsi + 1
450 :
451 31471 : if (use_full_chipsi) then
452 :
453 : buf_chipsi(ibuf_chipsi) = sum(cmplx(cg(1,icgb+1:icgb+npw),cg(2,icgb+1:icgb+npw),kind=dp)* &
454 : & paw_dmft%dpro(1:npw,iatom,ik)*paw_dmft%ylm(1:npw,im,lpawu+1,ik)* &
455 0 : & paw_dmft%bessel_int(1:npw,itypat,ik))
456 :
457 0 : if (prt_wan) then
458 :
459 0 : do ir=1,siz_wan
460 :
461 : ! Compute <Ylm|Psi_tilde>(r) = sum_g c_g * <Ylm|exp(j*(k+G)*(r+Rat))> / sqrt(ucvol)
462 : ! using exp(j*(k+G)*r) = 4*pi*sum_{lm} j**l * jl(|k+G|*r) * ylm(k+G) * ylm(theta,phi)
463 : ! (spherical harmonics expansion of planewave)
464 :
465 : paw_dmft%buf_psi(ibuf_psi+ir) = sum(paw_dmft%dpro(1:npw,iatom,ik) * paw_dmft%bessel(1:npw,ir,itypat,ik) * &
466 : & cmplx(cg(1,icgb+1:icgb+npw),cg(2,icgb+1:icgb+npw),kind=dp) * &
467 0 : & paw_dmft%ylm(1:npw,im,lpawu+1,ik))
468 :
469 : end do ! ir
470 :
471 : end if ! prt_wan
472 :
473 0 : if (verif_band) then
474 0 : do iproj=1,nproju
475 0 : buf_chipsi(ibuf_chipsi) = buf_chipsi(ibuf_chipsi) + cwprj(iproj,im)*paw_dmft%phimtphi_int(iproj,itypat)
476 0 : if (prt_wan) paw_dmft%buf_psi(ibuf_psi+1:ibuf_psi+siz_paw) = paw_dmft%buf_psi(ibuf_psi+1:ibuf_psi+siz_paw) + &
477 0 : & cwprj(iproj,im)*paw_dmft%phimtphi(1:siz_paw,iproj,itypat)
478 : end do ! iproj
479 : end if ! verif
480 :
481 : else
482 :
483 : ! In that case, simply assume |Psi> = \sum_{proja} <P_a|Psi><Chi|phi_a> (only true inside the PAW sphere
484 : ! and if your PAW basis is complete)
485 : ! Do not use DOT_PRODUCT
486 75893 : if (verif_band) buf_chipsi(ibuf_chipsi) = buf_chipsi(ibuf_chipsi) + sum(cwprj(1:nproju,im)*paw_dmft%phi_int(1:nproju,itypat))
487 :
488 : end if ! use_full_chipsi
489 :
490 51146 : ibuf_psi = ibuf_psi + siz_wan
491 :
492 : end do ! im
493 : end do ! iatom
494 12214 : icgb = icgb + npw
495 : end do ! ispinor
496 : end do ! iband
497 :
498 633 : icprj = icprj + nband_k_cprj*nspinor
499 993 : icg = icg + nband_k*npw*nspinor
500 : end do ! ikpt
501 : end do ! isppol
502 :
503 : !do isppol=1,nsppol
504 : !do ikpt=1,nkpt
505 : !do ispinor=1,my_nspinor
506 : !write(std_out,*) "psichi integers",isppol,ikpt,ispinor
507 : !write(std_out,*) "psichi IB3 iAT1 IM1",&
508 : !& real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1))
509 : !
510 : !enddo
511 : !enddo
512 : !enddo
513 : !call abi_abort('COLL')
514 : !if (abs(pawprtvol) >= 3) then
515 : ! write(message,*) "chinorm used here =",chinorm
516 : ! call wrtout(std_out,message,'COLL')
517 : !end if
518 :
519 : !deallocate temporary cwaveprj/cprj storage
520 97 : call pawcprj_free(cwaveprj)
521 423 : ABI_FREE(cwaveprj)
522 :
523 97 : ABI_FREE(cwprj)
524 97 : ABI_FREE(chipsi_tmp)
525 :
526 : !==========================================================================
527 : !********************* Gather information for MPI before printing
528 : !==========================================================================
529 :
530 : !call xmpi_barrier(comm_kpt)
531 97 : if (paral_kgb == 1 .and. nproc_band > 1) then
532 : ! Build sum over band processors
533 40 : call xmpi_sum(buf_chipsi(:),comm_band,ierr)
534 : end if
535 97 : call xmpi_allgatherv(buf_chipsi(:),siz_buf,buf_chipsi_tot(:),recvcounts(:),displs(:),comm_kpt,ierr)
536 :
537 97 : ABI_FREE(displs)
538 97 : ABI_FREE(recvcounts)
539 97 : ABI_FREE(buf_chipsi)
540 :
541 : ! Reorder the chipsi since the kpts can have any arbitrary distribution over the different MPI processes
542 97 : ibuf_chipsi = 0
543 290 : do irank=0,nproc_spkpt-1
544 581 : do isppol=1,nsppol
545 3562 : do ikpt=1,nkpt
546 :
547 3078 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
548 3078 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,irank)) cycle
549 :
550 11661 : do ibandc=1,mbandc
551 23638 : do ispinor=1,nspinor
552 47216 : do iatom=1,natom
553 26656 : lpawu = paw_dmft%lpawu(iatom)
554 26656 : if (lpawu == -1) cycle
555 11960 : ndim = 2*lpawu + 1
556 70536 : paw_dmft%chipsi(1+(ispinor-1)*ndim:ispinor*ndim,ibandc,ikpt,isppol,iatom) = buf_chipsi_tot(ibuf_chipsi+1:ibuf_chipsi+ndim)
557 25008 : ibuf_chipsi = ibuf_chipsi + ndim
558 : end do ! iatom
559 : end do ! ispinor
560 : end do ! ibandc
561 :
562 : end do ! ikpt
563 : end do ! isppol
564 : end do ! irank
565 :
566 : !call xmpi_barrier(comm_kpt)
567 :
568 97 : ABI_FREE(buf_chipsi_tot)
569 :
570 : !do isppol=1,nsppol
571 : !do ikpt=1,nkpt
572 : !do ibandc=1,paw_dmft%mbandc
573 : !do ispinor=1,my_nspinor
574 : !write(std_out,*) "psichigather",isppol,ikpt,ibandc,&
575 : !& real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2+&
576 : !& imag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2
577 : !
578 : !enddo
579 : !enddo
580 : !enddo
581 : !enddo
582 :
583 : !if(mpi_enreg%me.eq.0) write(177,*) "end",psichi
584 : !if(mpi_enreg%me.eq.1) write(178,*) "end",psichi
585 : !if(mpi_enreg%me.eq.2) write(179,*) "end",psichi
586 :
587 : !==========================================================================
588 : !********* WRITE unnormalized chipsi in file for reference
589 : !==========================================================================
590 97 : if (paw_dmft%myproc == 0) then
591 43 : call chipsi_print(paw_dmft,pawtab(:))
592 : end if ! proc=0
593 :
594 : !********************* Check normalization and occupations ***************
595 : ! Only if the complete BZ is sampled (ie paw_dmft%kspectralfunc=0)
596 : !==========================================================================
597 97 : if (paw_dmft%dmft_kspectralfunc == 0) then
598 :
599 93 : call init_oper(paw_dmft,loc_norm_check,opt_ksloc=2)
600 93 : call chipsi_check(paw_dmft,dft_occup,loc_norm_check)
601 : !==========================================================================
602 : !*************** write checks *******************************************
603 : !==========================================================================
604 : !if (abs(pawprtvol) >= 3) then
605 : ! write(message,*) "normalization computed"
606 : ! call wrtout(std_out,message,'COLL')
607 : !end if
608 :
609 93 : write(message,'(2a)') ch10," == The DMFT orbitals are now projected on the correlated bands"
610 93 : call wrtout(std_out,message,'COLL')
611 :
612 93 : write(message,'(2a,i4)') ch10," == Check: Downfolded Occupations and Norm of unnormalized projected orbitals"
613 93 : call wrtout(std_out,message,'COLL')
614 :
615 93 : if (paw_dmft%dmftcheck >= 1) then
616 : ! print occupations
617 1 : write(message,'(2a,i4)') ch10,' ------ Unsymmetrized Occupations'
618 1 : call wrtout(std_out,message,'COLL')
619 :
620 1 : call print_matlu(dft_occup%matlu(:),natom,pawprtvol)
621 :
622 : ! print norms
623 1 : write(message,'(2a,i4)') ch10,' ------ Unsymmetrized Norm'
624 1 : call wrtout(std_out,message,'COLL')
625 :
626 1 : call print_matlu(loc_norm_check%matlu(:),natom,pawprtvol)
627 : end if ! dmftcheck>=1
628 :
629 : ! symmetrize and print occupations
630 93 : call sym_matlu(dft_occup%matlu(:),paw_dmft)
631 :
632 93 : write(message,'(2a,i4)') ch10,' ------ Symmetrized Occupations'
633 93 : call wrtout(std_out,message,'COLL')
634 :
635 93 : call print_matlu(dft_occup%matlu(:),natom,pawprtvol)
636 :
637 : ! symmetrize and print norms
638 93 : call sym_matlu(loc_norm_check%matlu(:),paw_dmft)
639 :
640 93 : write(message,'(2a,i4)') ch10,' ------ Symmetrized Norm'
641 93 : call wrtout(std_out,message,'COLL')
642 :
643 93 : call print_matlu(loc_norm_check%matlu(:),natom,pawprtvol)
644 :
645 : ! Tests density matrix DFT+U and density matrix computed here.
646 93 : if (paw_dmft%dmftcheck == 2 .or. paw_dmft%dmftbandi == 1) then
647 28 : ABI_MALLOC(matlu_temp,(natom))
648 6 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu(:),matlu_temp(:))
649 6 : isppol = 1
650 6 : ispinor = 1
651 6 : ispinor1 = 1
652 6 : nsploop = max(nsppol,nspinor**2)
653 16 : do iatom=1,natom
654 10 : lpawu = paw_dmft%lpawu(iatom)
655 10 : if (lpawu == -1) cycle
656 6 : ndim = 2*lpawu + 1
657 24 : do idijeff=1,nsploop
658 12 : if (nsploop <= 2) then
659 12 : isppol = idijeff
660 0 : else if (nsploop == 4) then
661 0 : ispinor = spinor_idxs(1,idijeff)
662 0 : ispinor1 = spinor_idxs(2,idijeff)
663 : else
664 0 : ABI_BUG(" BUG in datafordmft: nsploop should be equal to 2 or 4")
665 : end if ! nsploop
666 12 : if (my_nspinor == 2) matlu_temp(iatom)%mat(1+(ispinor-1)*ndim:ispinor*ndim,1+(ispinor1-1)*ndim:ispinor1*ndim,isppol) = &
667 0 : & cmplx(paw_ij(iatom)%noccmmp(1,1:ndim,1:ndim,idijeff),paw_ij(iatom)%noccmmp(2,1:ndim,1:ndim,idijeff),kind=dp)
668 12 : if (my_nspinor == 1) matlu_temp(iatom)%mat(1+(ispinor-1)*ndim:ispinor*ndim,1+(ispinor1-1)*ndim:ispinor1*ndim,isppol) = &
669 486 : & cmplx(paw_ij(iatom)%noccmmp(1,1:ndim,1:ndim,idijeff),zero,kind=dp)
670 : end do ! idijeff
671 : end do ! iatom
672 6 : if (paw_dmft%dmftcheck == 2) option = 1
673 6 : if (paw_dmft%dmftcheck <= 1) option = 0
674 : call diff_matlu("DFT+U density matrix from INPUT wfk",&
675 : & "Direct calculation of density matrix with chipsi from DIAGONALIZED wfk",&
676 6 : & matlu_temp(:),dft_occup%matlu(:),natom,option,tol3,ierr) !tol1 tol2 tol3
677 6 : if (ierr == -1) then
678 2 : write(message,'(10a)') ch10,&
679 2 : & ' -> These two quantities should agree if three conditions are fullfilled',ch10,&
680 2 : & ' - input wavefunctions come from the same Hamiltonian (e.g LDA/GGA)',ch10,&
681 2 : & ' - dmatpuopt is equal to 1',ch10,&
682 2 : & ' - all valence states are in the valence',ch10,&
683 4 : & ' (for experts users: it is not compulsory that these conditions are fullfilled)'
684 2 : call wrtout(std_out,message,'COLL')
685 : end if
686 : ! write(message,'(2a)') ch10,&
687 : ! & ' ***** => Calculations of density matrices with projections and in DFT+U are coherent****'
688 : ! call wrtout(std_out,message,'COLL')
689 :
690 6 : call destroy_matlu(matlu_temp(:),natom)
691 16 : ABI_FREE(matlu_temp)
692 : else
693 87 : write(message,'(2a)') ch10,&
694 174 : & ' Warning: Consistency of density matrices computed from projection has not been checked: use dmftcheck>=2 '
695 87 : call wrtout(std_out,message,'COLL')
696 : end if
697 :
698 93 : call destroy_oper(loc_norm_check)
699 : end if ! dmft_kspectralfunc=0
700 :
701 194 : if (present(nbandkss)) then
702 97 : if ((me_kpt == 0 .and. nbandkss /= 0) .or. (paw_dmft%dmft_kspectralfunc == 1)) then
703 : ! opt_renorm=1 ! if ucrpa==1, no need for individual orthonormalization
704 4 : opt_renorm = merge(2,3,dtset%ucrpa>=1.or.paw_dmft%dmft_kspectralfunc==1)
705 4 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) opt_renorm = paw_dmft%dmft_wanorthnorm
706 4 : call chipsi_renormalization(paw_dmft,opt=opt_renorm)
707 4 : if (paw_dmft%myproc == 0) then
708 1 : call chipsi_print(paw_dmft,pawtab(:))
709 : end if
710 : end if ! proc=me
711 : end if
712 :
713 : CONTAINS
714 :
715 : !!***
716 : !!****f* m_datafordmft/chipsi_check
717 : !! NAME
718 : !! chipsi_check
719 : !!
720 : !! FUNCTION
721 : !! Check chipsi: compute norm and occupations
722 : !!
723 : !! INPUTS
724 : !! paw_dmft <type(paw_dmft)>=paw data for the self-consistency
725 : !!
726 : !! OUTPUTS:
727 : !! xocc_check: density matrix
728 : !! xnorm_check: matrix of norms
729 : !!
730 : !! SIDE EFFECTS
731 : !!
732 : !! SOURCE
733 :
734 93 : subroutine chipsi_check(paw_dmft,xocc_check,xnorm_check)
735 :
736 : !Arguments ------------------------------------
737 : type(paw_dmft_type), intent(in) :: paw_dmft
738 : type(oper_type), intent(inout) :: xnorm_check,xocc_check
739 : !Local variables ------------------------------------
740 93 : real(dp), allocatable :: occ_dft(:,:,:)
741 : ! *********************************************************************
742 :
743 465 : ABI_MALLOC(occ_dft,(mbandc,nkpt,nsppol))
744 :
745 93 : band_index = 0
746 224 : do isppol=1,nsppol
747 1302 : do ikpt=1,nkpt
748 1078 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
749 1078 : ibandc = 0
750 26118 : do iband=1,nband_k
751 26118 : if (paw_dmft%band_in(iband)) then
752 9852 : ibandc = ibandc + 1
753 9852 : occ_dft(ibandc,ikpt,isppol) = occ(iband+band_index)
754 : end if
755 : end do ! iband
756 1209 : band_index = band_index + nband_k
757 : end do ! ikpt
758 : end do ! isppol
759 :
760 4941 : if (nsppol == 1 .and. my_nspinor == 1) occ_dft(:,:,:) = occ_dft(:,:,:) * half
761 :
762 93 : call downfold_oper(xnorm_check,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=2)
763 93 : call xmpi_matlu(xnorm_check%matlu(:),natom,paw_dmft%distrib%comm_kpt)
764 : call downfold_oper(xocc_check,paw_dmft,procb=paw_dmft%distrib%procb(:), &
765 93 : & iproc=paw_dmft%distrib%me_kpt,option=3,op_ks_diag=occ_dft(:,:,:))
766 93 : call xmpi_matlu(xocc_check%matlu(:),natom,paw_dmft%distrib%comm_kpt)
767 :
768 93 : ABI_FREE(occ_dft)
769 :
770 93 : end subroutine chipsi_check
771 : !DBG_EXIT("COLL")
772 : !!***
773 : end subroutine datafordmft
774 : !!***
775 :
776 : !!****f* m_datafordmft/chipsi_print
777 : !! NAME
778 : !! chipsi_print
779 : !!
780 : !! FUNCTION
781 : !! Print chipsi for reference
782 : !!
783 : !! INPUTS
784 : !! paw_dmft <type(paw_dmft)>=paw data for the self-consistency
785 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
786 : !!
787 : !! SIDE EFFECTS
788 : !! print chipsi in forlb.ovlp
789 : !!
790 : !! SOURCE
791 :
792 44 : subroutine chipsi_print(paw_dmft,pawtab)
793 :
794 : !Arguments ------------------------------------
795 : type(paw_dmft_type), intent(in) :: paw_dmft
796 : type(pawtab_type), intent(in) :: pawtab(paw_dmft%ntypat)
797 : !Local variables ------------------------------------
798 : integer :: iatom,iband,ibandc,ikpt,im,im1,ispinor,isppol
799 : integer :: itypat,lpawu,nband_k,ndim,unt
800 : logical :: t2g,x2my2d
801 : character(len=500) :: msg
802 : integer, parameter :: mt2g(3) = (/1,2,4/)
803 : ! *********************************************************************
804 :
805 44 : t2g = (paw_dmft%dmft_t2g == 1)
806 44 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
807 :
808 44 : if (open_file('forlb.ovlp',msg,newunit=unt,form='formatted',status='unknown') /= 0) ABI_ERROR(msg)
809 44 : rewind(unt)
810 :
811 : ! Header for calc_uCRPA.F90
812 44 : write(unt,*) "# isppol nspinor natom m Re(<chi|psi>) Im(<chi|psi>)"
813 132 : if (count(pawtab(:)%lpawu /= -1) == 1) then
814 132 : do itypat=1,paw_dmft%ntypat
815 88 : lpawu = pawtab(itypat)%lpawu
816 88 : if (lpawu == -1) cycle
817 88 : if (t2g) then
818 17 : write(unt,*) "l= ",1,itypat
819 27 : else if (x2my2d) then
820 0 : write(unt,*) "l= ",0,itypat
821 : else
822 27 : write(unt,*) "l= ",lpawu,itypat
823 : end if
824 : end do ! itypat
825 : else
826 0 : write(unt,*) "More than one correlated species"
827 : end if
828 :
829 44 : write(unt,*) "Bands ",paw_dmft%dmftbandi,paw_dmft%dmftbandf
830 :
831 111 : do isppol=1,paw_dmft%nsppol
832 605 : do ikpt=1,paw_dmft%nkpt
833 : ! rewind(1023)
834 494 : write(unt,'(a6,2x,i6)') "ikpt =",ikpt
835 494 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*paw_dmft%nkpt)
836 494 : ibandc = 0
837 11797 : do iband=1,nband_k
838 11236 : if (paw_dmft%band_in(iband)) then
839 4374 : ibandc = ibandc + 1
840 4374 : write(unt,'(a8,2x,i6)') " iband =",iband
841 : else
842 : cycle
843 : end if
844 9306 : do ispinor=1,paw_dmft%nspinor
845 26448 : do iatom=1,paw_dmft%natom
846 10774 : lpawu = paw_dmft%lpawu(iatom)
847 10774 : if (lpawu == -1) cycle
848 5150 : ndim = 2*lpawu + 1
849 36154 : do im=1,ndim
850 26566 : if (t2g) then
851 3060 : im1 = mt2g(im)
852 23506 : else if (x2my2d) then
853 0 : im1 = 5
854 : else
855 23506 : im1 = im
856 : end if
857 :
858 26566 : write(unt,'(4i6,3x,2f23.15)') isppol,ispinor,iatom,im1,&
859 26566 : & dble(paw_dmft%chipsi(im+(ispinor-1)*ndim,ibandc,ikpt,isppol,iatom)),&
860 63906 : & aimag(paw_dmft%chipsi(im+(ispinor-1)*ndim,ibandc,ikpt,isppol,iatom))
861 : end do !im
862 : end do ! iatom
863 : end do ! ispinor
864 : end do !iband
865 : end do !ikpt
866 : end do ! isppol
867 : ! write(unt,*) "Fermi level (in Ryd)="
868 : ! write(unt,*) fermie*two
869 44 : close(unt)
870 :
871 44 : end subroutine chipsi_print
872 : !!***
873 :
874 : !!****f* m_datafordmft/compute_levels
875 : !! NAME
876 : !! compute_levels
877 : !!
878 : !! FUNCTION
879 : !! Compute correlated electronic levels for ctqmc
880 : !!
881 : !! INPUTS
882 : !! hdc= double counting
883 : !! paw_dmft <type(paw_dmft)>=paw data for the self-consistency
884 : !!
885 : !! OUTPUT
886 : !! energy_level= local electronic levels
887 : !! nondiag= true if the levels are not diagonal
888 : !!
889 : !! NOTES
890 : !!
891 : !! SOURCE
892 :
893 209 : subroutine compute_levels(energy_level,hdc,paw_dmft,nondiag)
894 :
895 : !Arguments ------------------------------------
896 : type(oper_type), intent(in) :: hdc
897 : type(paw_dmft_type), intent(in) :: paw_dmft
898 : type(oper_type), intent(inout) :: energy_level
899 : logical, optional, intent(out) :: nondiag
900 : !Local variables ------------------------------
901 : integer :: iatom,im,isppol,lpawu,natom,ndim,nspinor,nsppol
902 : character(len=13) :: tag
903 : character(len=500) :: message
904 : !type
905 209 : type(matlu_type), allocatable :: levels_temp(:),magnfield(:)
906 : !************************************************************************
907 :
908 209 : natom = paw_dmft%natom
909 209 : nspinor = paw_dmft%nspinor
910 209 : nsppol = paw_dmft%nsppol
911 209 : if (present(nondiag)) nondiag = .false.
912 :
913 : !======================================================================
914 : !Compute atomic levels from projection of \epsilon_{nks} and symmetrize
915 : !======================================================================
916 :
917 : call downfold_oper(energy_level,paw_dmft,option=3,procb=paw_dmft%distrib%procb(:), &
918 209 : & iproc=paw_dmft%distrib%me_kpt,op_ks_diag=paw_dmft%eigen_dft(:,:,:))
919 209 : call xmpi_matlu(energy_level%matlu(:),natom,paw_dmft%distrib%comm_kpt)
920 : ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels before sym and only DFT"
921 : ! call wrtout(std_out,message,'COLL')
922 : ! call print_matlu(energy_level%matlu,natom,1)
923 943 : do iatom=1,natom
924 734 : lpawu = paw_dmft%lpawu(iatom)
925 734 : if (lpawu == -1) cycle
926 257 : ndim = paw_dmft%nspinor * (2*lpawu+1)
927 659 : do isppol=1,paw_dmft%nsppol
928 2561 : do im=1,ndim
929 2304 : energy_level%matlu(iatom)%mat(im,im,isppol) = energy_level%matlu(iatom)%mat(im,im,isppol) - paw_dmft%fermie
930 : end do ! im
931 : end do ! isppol
932 13117 : energy_level%matlu(iatom)%mat(:,:,:) = energy_level%matlu(iatom)%mat(:,:,:) - hdc%matlu(iatom)%mat(:,:,:)
933 : ! write(std_out,*) "DC,fermie",hdc%matlu(iatom)%mat(1,1,1,1,1),paw_dmft%fermie
934 : end do ! iatom
935 :
936 209 : call sym_matlu(energy_level%matlu(:),paw_dmft)
937 209 : if (present(nondiag)) then
938 140 : call checkdiag_matlu(energy_level%matlu(:),natom,tol7,nondiag)
939 : end if
940 :
941 209 : write(tag,'(f13.5)') paw_dmft%fermie
942 209 : write(message,'(a,2x,2a)') ch10," == Print Energy levels in cubic basis for Fermi Level = ",adjustl(tag)
943 209 : call wrtout(std_out,message,'COLL')
944 : !call print_oper(energy_level,1,paw_dmft,1)
945 209 : call print_matlu(energy_level%matlu(:),natom,1)
946 :
947 : !========================================================
948 : !Add Zeeman contibutions to local energy levels in slm
949 : !========================================================
950 209 : if(paw_dmft%dmft_magnfield .eq. 2 .and. nspinor .eq. 1) then
951 32 : ABI_MALLOC(magnfield,(natom))
952 24 : ABI_MALLOC(levels_temp,(natom))
953 :
954 8 : write(message,'(a,2x,2a)') ch10, " == Adding Zeeman contribution to local energy levels"
955 8 : call wrtout(std_out,message,'COLL')
956 :
957 : !Apply-m.bz only in the local Hamiltonian
958 8 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,magnfield)
959 8 : call init_matlu(natom,nspinor,nsppol,paw_dmft%lpawu,levels_temp)
960 8 : call copy_matlu(energy_level%matlu,levels_temp,natom)
961 :
962 8 : call magnfield_matlu(magnfield,natom,paw_dmft%dmft_magnfield_b,1)
963 : !call print_matlu(magnfield,natom,1)
964 8 : call add_matlu(levels_temp,magnfield,energy_level%matlu,natom,-1)
965 8 : call print_matlu(energy_level%matlu,natom,1)
966 8 : call destroy_matlu(magnfield,natom)
967 8 : call destroy_matlu(levels_temp,natom)
968 :
969 16 : ABI_FREE(magnfield)
970 16 : ABI_FREE(levels_temp)
971 : endif
972 :
973 209 : end subroutine compute_levels
974 : !!***
975 :
976 : !!****f* m_datafordmft/chipsi_renormalization
977 : !! NAME
978 : !! chipsi_renormalization
979 : !!
980 : !! FUNCTION
981 : !! Orthonormalize chipsi.
982 : !!
983 : !! INPUTS
984 : !! paw_dmft = data for DFT+DMFT calculations.
985 : !! opt = 2 : orthonormalize all the atoms at the same time,
986 : !! and for each individual kpt
987 : !! = 3 (default) : orthonormalize the sum over all kpt,
988 : !! and for each individual atom
989 : !!
990 : !! OUTPUT
991 : !! paw_dmft%chipsi((2*maxlpawu+1)*nspinor,mbandc,nkpt,nsppol,natom):
992 : !! orthonormalized projections <Chi|Psi>
993 : !!
994 : !! NOTES
995 : !!
996 : !! SOURCE
997 :
998 96 : subroutine chipsi_renormalization(paw_dmft,opt)
999 :
1000 : !Arguments ------------------------------------
1001 : type(paw_dmft_type), intent(inout) :: paw_dmft
1002 : integer, optional, intent(in) :: opt
1003 : !Local variables ------------------------------
1004 : integer :: jkpt,natom,nkpt,option
1005 : real(dp) :: pawprtvol
1006 96 : type(oper_type) :: norm,oper_temp
1007 : character(len=500) :: message
1008 96 : real(dp), allocatable :: wtk_tmp(:)
1009 : character(len=5) :: tag
1010 : ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:)
1011 : !************************************************************************
1012 :
1013 : DBG_ENTER("COLL")
1014 :
1015 96 : option = 3
1016 96 : if (present(opt)) then
1017 96 : if (opt == 2 .or. opt == 3) option = opt
1018 : end if
1019 96 : pawprtvol = 2
1020 :
1021 96 : natom = paw_dmft%natom
1022 96 : nkpt = paw_dmft%nkpt
1023 :
1024 : !== Normalize psichi
1025 : !if (option == 1) then
1026 : ! ====================================
1027 : ! == simply renormalize psichi =======
1028 : ! ====================================
1029 : ! write(message,'(2a)') ch10," Psichi are renormalized "
1030 : ! call wrtout(std_out, message,'COLL')
1031 : ! do isppol=1,nsppol
1032 : ! do ikpt=1,nkpt
1033 : ! do ib=1,mbandc
1034 : ! do iatom=1,natom
1035 : ! if(paw_dmft%lpawu(iatom).ne.-1) then
1036 : ! ndim=2*paw_dmft%lpawu(iatom)+1
1037 : ! do im=1,ndim
1038 : ! do ispinor=1,nspinor
1039 : ! write(std_out,*) "psichi1",paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)
1040 : ! paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)= &
1041 : !& paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)/ &
1042 : !& sqrt(real(norm%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)))
1043 : ! end do ! ispinor
1044 : ! end do ! im
1045 : ! end if
1046 : ! end do ! iatom
1047 : ! end do ! ib
1048 : ! end do ! ikpt
1049 : ! end do ! isppol
1050 : ! todo_ab introduce correct orthonormalization in the general case.
1051 :
1052 96 : write(message,'(6a)') ch10, &
1053 96 : & ' =================================================== ',&
1054 96 : & ch10,' == The DMFT orbitals will now be orthonormalized == ',&
1055 192 : & ch10,' =================================================== '
1056 96 : call wrtout(std_out,message,'COLL')
1057 :
1058 96 : if (option == 2) then ! option==2
1059 : ! ====================================
1060 : ! == renormalize k-point after k-point
1061 : ! ====================================
1062 :
1063 17 : write(tag,'(i5)') nkpt
1064 17 : write(message,'(3a)') ch10,' Number of k-points: ',adjustl(tag)
1065 17 : call wrtout(std_out,message,'COLL')
1066 321 : do jkpt=1,nkpt ! jkpt
1067 304 : write(tag,'(i5)') jkpt
1068 304 : write(message,'(3a)') ch10,' == Orthonormalization for k-point: ',adjustl(tag)
1069 304 : call wrtout(std_out,message,'COLL')
1070 304 : if (paw_dmft%distrib%procb(jkpt) /= paw_dmft%distrib%me_kpt) cycle
1071 321 : call normalizechipsi(1,paw_dmft,jkpt=jkpt)
1072 : end do ! jkpt
1073 17 : write(message,'(2a)') ch10,' ===== Every k-point has been orthonormalized'
1074 17 : call wrtout(std_out,message,'COLL')
1075 :
1076 : else if (option == 3) then ! option==3
1077 : ! ====================================
1078 : ! == renormalize the sum over k-points
1079 : ! ====================================
1080 79 : write(message,'(6a)') ch10, &
1081 79 : & ' ====================================== ',&
1082 79 : & ch10,' == Orthonormalization for all k-points == ',&
1083 158 : & ch10,' ====================================== '
1084 79 : call wrtout(std_out,message,'COLL')
1085 79 : call normalizechipsi(nkpt,paw_dmft)
1086 :
1087 : end if ! option
1088 :
1089 : ! Gather contribution from each CPU
1090 96 : call chipsi_gather(paw_dmft)
1091 :
1092 : !== Change back repr for norm
1093 :
1094 : !===============================================
1095 : !== Compute norm with new chipsi
1096 : !===============================================
1097 :
1098 96 : write(message,'(2a)') ch10,' ===== Compute new norm after renormalization'
1099 96 : call wrtout(std_out,message,'COLL')
1100 96 : call init_oper(paw_dmft,oper_temp,opt_ksloc=2)
1101 96 : call identity_oper(oper_temp,2)
1102 :
1103 96 : if (paw_dmft%dmft_kspectralfunc == 1 .and. &
1104 : & ((paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7) .or. paw_dmft%kptopt < 0)) then
1105 12 : ABI_MALLOC(wtk_tmp,(nkpt))
1106 48 : wtk_tmp(:) = one
1107 4 : call init_oper(paw_dmft,norm,nkpt=1,wtk=wtk_tmp(:),opt_ksloc=2)
1108 48 : do jkpt=1,nkpt ! jkpt
1109 44 : norm%shiftk = jkpt - 1
1110 44 : call downfold_oper(norm,paw_dmft,option=2)
1111 44 : write(tag,'(i5)') jkpt
1112 : write(message,'(3a)') &
1113 44 : & ch10," == Check: Overlap after orthonormalization for k-point ",adjustl(tag)
1114 44 : call wrtout(std_out,message,'COLL')
1115 44 : call print_matlu(norm%matlu(:),natom,prtopt=1)
1116 : !== Check that norm is now the identity
1117 : call diff_matlu('Overlap after orthonormalization','Identity',&
1118 48 : & norm%matlu(:),oper_temp%matlu(:),natom,1,tol6,zero_or_one=1)
1119 : end do ! jkpt
1120 4 : ABI_FREE(wtk_tmp)
1121 : else !dmft_kspectralfunc
1122 92 : write(message,'(2a)') ch10,' ===== Starting downfold'
1123 92 : call wrtout(std_out,message,'COLL')
1124 92 : call init_oper(paw_dmft,norm,opt_ksloc=2)
1125 92 : call downfold_oper(norm,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=2)
1126 92 : call xmpi_matlu(norm%matlu(:),natom,paw_dmft%distrib%comm_kpt)
1127 92 : write(message,'(2a)') ch10,' ===== Finished downfold'
1128 92 : call wrtout(std_out,message,'COLL')
1129 :
1130 : !== Print unsymmetrized norm%matlu with new chipsi
1131 : if (pawprtvol > 2) then
1132 : write(message,'(4a,2a)') &
1133 : & ch10," == Check: Overlap with renormalized chipsi without symmetrization is == "
1134 : call wrtout(std_out,message,'COLL')
1135 : call print_matlu(norm%matlu(:),natom,prtopt=1)
1136 : end if ! pawprtvol>2
1137 :
1138 : !== Symmetrize norm%matlu with new chipsi
1139 92 : call sym_matlu(norm%matlu(:),paw_dmft)
1140 :
1141 : !== Print symmetrized norm%matlu with new chipsi
1142 : if (pawprtvol > 2) then
1143 : write(message,'(4a,2a)') &
1144 : & ch10," == Check: Overlap with renormalized chipsi and symmetrization is =="
1145 : call wrtout(std_out,message,'COLL')
1146 : call print_matlu(norm%matlu(:),natom,prtopt=1,opt_diag=-1)
1147 : end if
1148 :
1149 : !== Check that norm is now the identity
1150 92 : call diff_oper('Overlap after renormalization','Identity',norm,oper_temp,0,tol6)
1151 :
1152 : end if ! dmft_kspectralfunc
1153 :
1154 96 : call destroy_oper(norm)
1155 96 : call destroy_oper(oper_temp)
1156 :
1157 96 : paw_dmft%lchipsiortho = 1
1158 :
1159 : DBG_EXIT("COLL")
1160 :
1161 : CONTAINS
1162 : !===========================================================
1163 : !!***
1164 :
1165 : !!****f* chipsi_renormalization/normalizechipsi
1166 : !! NAME
1167 : !! normalizechipsi
1168 : !!
1169 : !! FUNCTION
1170 : !! Orthonormalize chipsi
1171 : !!
1172 : !! INPUTS
1173 : !! nkpt = number of kpt
1174 : !! paw_dmft = data for DFT+DMFT calculations.
1175 : !! jkpt = if present, index of the kpt to be orthonormalized
1176 : !!
1177 : !! SIDE EFFECTS
1178 : !!
1179 : !! SOURCE
1180 :
1181 170 : subroutine normalizechipsi(nkpt,paw_dmft,jkpt)
1182 :
1183 : !Arguments ------------------------------------
1184 : integer, intent(in) :: nkpt
1185 : integer, optional, intent(in) :: jkpt
1186 : type(paw_dmft_type), intent(inout) :: paw_dmft
1187 : !Local variables ------------------------------
1188 : integer :: dimoverlap,dum,iatom,ib,ikpt,im,im1,isppol,itot,itot1,lpawu,mbandc
1189 : integer :: natom,ndim,ndim_max,nspinor,nsppol,pawprtvol,unt
1190 : logical :: lexist
1191 170 : type(oper_type) :: norm1,norm2,norm3
1192 170 : real(dp), allocatable :: si(:),sr(:)
1193 170 : complex(dp), allocatable :: chipsivect(:,:),largeoverlap(:,:),mat_tmp(:,:)
1194 : character(len=1) :: tag_is
1195 : character(len=4) :: tag_at
1196 : character(len=5) :: tag
1197 : character(len=500) :: message
1198 : character(len=fnlen) :: tmpfil
1199 : ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:)
1200 : !complex(dp), allocatable :: wan(:,:,:),sqrtmatinv(:,:),wanall(:)
1201 : !type(coeff2c_type), allocatable :: overlap(:)
1202 : !************************************************************************
1203 :
1204 170 : mbandc = paw_dmft%mbandc
1205 170 : natom = paw_dmft%natom
1206 170 : nspinor = paw_dmft%nspinor
1207 170 : nsppol = paw_dmft%nsppol
1208 170 : ndim_max = (2*paw_dmft%maxlpawu+1) * nspinor
1209 170 : pawprtvol = 3
1210 :
1211 170 : if (nkpt /= 1 .and. present(jkpt)) ABI_BUG('BUG in chipsi_normalization')
1212 :
1213 : ! iortho=1
1214 : ! write(6,*) "nkpt, iortho",nkpt,iortho
1215 : !if (natomcor>1) iortho=2
1216 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1217 : ! First case: usual case (this numerically guarantees downfold(upfold)=Id only for one atom and nkpt=1)
1218 98 : if (.not. present(jkpt)) then ! .and.iortho==1) then
1219 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1220 : ! *********************************************************************
1221 :
1222 79 : call init_oper(paw_dmft,norm1,opt_ksloc=2)
1223 :
1224 79 : if (paw_dmft%dmft_kspectralfunc == 1 .and. paw_dmft%kptopt < 0) then ! Read overlap from file
1225 0 : if (paw_dmft%myproc == 0) then
1226 0 : ABI_MALLOC(sr,(ndim_max))
1227 0 : ABI_MALLOC(si,(ndim_max))
1228 0 : do iatom=1,natom
1229 0 : lpawu = paw_dmft%lpawu(iatom)
1230 0 : if (lpawu == -1) cycle
1231 0 : ndim = nspinor * (2*lpawu+1)
1232 0 : call int2char4(iatom,tag_at)
1233 0 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
1234 0 : do isppol=1,nsppol
1235 0 : write(tag_is,'(i1)') isppol
1236 0 : tmpfil = trim(paw_dmft%filnamei)//"_DMFTOVERLAP_iatom"//tag_at//"_isppol"//tag_is
1237 0 : inquire(file=trim(tmpfil),exist=lexist)
1238 0 : if (.not. lexist) ABI_ERROR("File "//trim(tmpfil)//" does not exist !")
1239 0 : if (open_file(trim(tmpfil),message,newunit=unt) /= 0) ABI_ERROR(message)
1240 0 : do im=1,ndim
1241 0 : read(unt,'(60(2x,e25.17e3))') (sr(im1),si(im1),im1=1,ndim)
1242 0 : norm1%matlu(iatom)%mat(im,:,isppol) = cmplx(sr(1:ndim),si(1:ndim),kind=dp)
1243 : end do ! im
1244 0 : close(unt)
1245 : end do ! isppol
1246 : end do ! iatom
1247 0 : ABI_FREE(sr)
1248 0 : ABI_FREE(si)
1249 : end if ! myproc=0
1250 0 : call xmpi_matlu(norm1%matlu(:),natom,paw_dmft%spacecomm,0,2)
1251 : else ! Compute overlap
1252 79 : call downfold_oper(norm1,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=2)
1253 79 : call xmpi_matlu(norm1%matlu(:),natom,paw_dmft%distrib%comm_kpt)
1254 79 : if (nkpt > 1) then
1255 72 : call sym_matlu(norm1%matlu(:),paw_dmft)
1256 : end if
1257 : end if ! kspectralfunc
1258 :
1259 : if (pawprtvol > 2) then
1260 79 : write(message,'(2a)') ch10,' - Print current norm and overlap (before orthonormalization)'
1261 79 : call wrtout(std_out,message,'COLL')
1262 79 : call print_matlu(norm1%matlu(:),natom,prtopt=1,opt_exp=1)
1263 : end if
1264 :
1265 : ! build large overlap matrix
1266 : !write(message,'(2a)') ch10,' - Overlap (before orthonormalization) -'
1267 : !call wrtout(std_out,message,'COLL')
1268 :
1269 : ! ==-------------------------------------
1270 : ! == Start loop over atoms
1271 :
1272 : !ABI_MALLOC(overlap,(natom))
1273 : !do iatom=1,natom
1274 : ! if(paw_dmft%lpawu(iatom).ne.-1) then
1275 : ! ndim=2*paw_dmft%lpawu(iatom)+1
1276 : ! tndim=nsppol*nspinor*ndim
1277 : ! ABI_MALLOC(overlap(iatom)%value,(tndim,tndim))
1278 : ! overlap(iatom)%value=czero
1279 : ! end if
1280 : !end do
1281 : ! ==-------------------------------------
1282 :
1283 : ! built large overlap matrix
1284 : !write(message,'(2a)') ch10,' - Overlap (before orthonormalization) -'
1285 : !call wrtout(std_out,message,'COLL')
1286 : !call gather_matlu(norm1%matlu,overlap,cryst_struc%natom,option=1,prtopt=1)
1287 : !call destroy_oper(norm1)
1288 :
1289 385 : do iatom=1,natom
1290 306 : lpawu = paw_dmft%lpawu(iatom)
1291 306 : if (lpawu == -1) cycle
1292 91 : ndim = nspinor * (2*lpawu+1)
1293 364 : ABI_MALLOC(mat_tmp,(ndim,mbandc))
1294 :
1295 : ! == Compute Inverse Square root of overlap : O^{-0.5}
1296 : !do im=1,tndim
1297 : ! do im1=1,tndim
1298 : ! !write(message,'(a,1x,a,e21.14,a,e21.14,a)') "overlap", &
1299 : ! !"(",real(overlap(1)%value(im,im1)),",",aimag(overlap(1)%value(im,im1)),")"
1300 : ! write(6,*) "overlap",overlap(iatom)%value(im,im1)
1301 : ! enddo
1302 : !enddo
1303 : !stop
1304 : !call wrtout(std_out,message,'COLL')
1305 : !if(diag==0) then
1306 : !call invsqrt_matrix(overlap(iatom)%value,tndim,dum)
1307 : !sqrtmatinv=overlap(iatom)%value
1308 218 : do isppol=1,nsppol
1309 :
1310 : if ((paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) .and. paw_dmft%dmft_prt_maxent == 1 .and. paw_dmft%myproc == 0 &
1311 127 : & .and. paw_dmft%dmft_kspectralfunc == 0) then
1312 0 : call int2char4(iatom,tag_at)
1313 0 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
1314 0 : write(tag_is,'(i1)') isppol
1315 0 : if (open_file(trim(paw_dmft%filapp)//"_DMFTOVERLAP_iatom"//tag_at//"_isppol"//tag_is,message,newunit=unt) /= 0) ABI_ERROR(message)
1316 0 : do im=1,ndim
1317 0 : write(unt,'(60(2x,e25.17e3))') (norm1%matlu(iatom)%mat(im,im1,isppol),im1=1,ndim)
1318 : end do ! im
1319 0 : close(unt)
1320 : end if ! triqs and dmft_prt_maxent
1321 :
1322 : ! if(diag==0) then
1323 127 : call invsqrt_matrix(norm1%matlu(iatom)%mat(:,:,isppol),ndim,dum)
1324 : !else
1325 : ! do im1=1,tndim
1326 : ! do im=1,tndim
1327 : ! if (im==im1) then
1328 : ! norm1%matlu(iatom)%mat(im,im1,isppol)=cone/sqrt(norm1%matlu(iatom)%mat(im,im,isppol)
1329 : ! else
1330 : ! norm1%matlu(iatom)%mat(im,im1,isppol)=czero
1331 : ! end if
1332 : ! end do
1333 : ! end do
1334 : !end if
1335 :
1336 : ! == Apply O^{-0.5} on chipsi
1337 838 : do ikpt=1,nkpt
1338 620 : if (paw_dmft%distrib%procb(ikpt) /= paw_dmft%distrib%me_kpt) cycle
1339 : call abi_xgemm("n","n",ndim,mbandc,ndim,cone,norm1%matlu(iatom)%mat(:,:,isppol),ndim, &
1340 338 : & paw_dmft%chipsi(:,:,ikpt,isppol,iatom),ndim_max,czero,mat_tmp(:,:),ndim)
1341 24989 : paw_dmft%chipsi(1:ndim,:,ikpt,isppol,iatom) = mat_tmp(:,:)
1342 : end do ! ikpt
1343 : end do ! isppol
1344 :
1345 :
1346 : ! ABI_MALLOC(wan,(nsppol,nspinor,ndim))
1347 : ! write(std_out,*) mbandc,nsppol,nspinor,ndim
1348 : ! write(std_out,*) paw_dmft%psichi(1,1,1,1,1,1)
1349 : ! do ikpt=1,nkpt
1350 : ! do ib=1,mbandc
1351 : ! if(present(jkpt)) then
1352 : ! ikpt1=jkpt
1353 : ! else
1354 : ! ikpt1=ikpt
1355 : ! end if
1356 : ! jc=0
1357 : ! wan=czero
1358 : ! do isppol=1,nsppol
1359 : ! do ispinor=1,nspinor
1360 : ! do im=1,ndim
1361 : ! write(std_out,*) "psichi", paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)
1362 : ! jc=jc+1
1363 : ! jc1=0
1364 : ! do isppol1=1,nsppol
1365 : ! do ispinor1=1,nspinor
1366 : ! do im1=1,ndim
1367 : ! jc1=jc1+1
1368 : ! wan(isppol,ispinor,im)= wan(isppol,ispinor,im) &
1369 : !& + paw_dmft%psichi(isppol1,ikpt1,ib,ispinor1,iatom,im1)*sqrtmatinv(jc,jc1)
1370 : ! end do ! ispinor1
1371 : ! end do ! isppol1
1372 : ! end do ! im1
1373 : ! end do ! im
1374 : ! end do ! ispinor
1375 : ! end do ! isppol
1376 : ! do isppol=1,nsppol
1377 : ! do ispinor=1,nspinor
1378 : ! do im=1,ndim
1379 : ! paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)=wan(isppol,ispinor,im)
1380 : ! write(std_out,*) "psichi2", paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)
1381 : ! end do ! ispinor
1382 : ! end do ! isppol
1383 : ! end do ! im
1384 : ! end do ! ib
1385 : ! end do ! ikpt
1386 : ! ABI_FREE(wan)
1387 : ! ABI_FREE(sqrtmatinv)
1388 :
1389 385 : ABI_FREE(mat_tmp)
1390 :
1391 : end do ! iatom
1392 :
1393 : ! == End loop over atoms
1394 : ! ==-------------------------------------
1395 :
1396 : ! ======================================================================
1397 : ! == Check norm with new chipsi.
1398 : ! ======================================================================
1399 :
1400 79 : call downfold_oper(norm1,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=2)
1401 79 : call xmpi_matlu(norm1%matlu(:),natom,paw_dmft%distrib%comm_kpt)
1402 :
1403 79 : if (nkpt > 1) then
1404 72 : call sym_matlu(norm1%matlu(:),paw_dmft)
1405 : end if
1406 :
1407 : if (pawprtvol > 2) then
1408 79 : write(message,'(2a)') ch10,' - Print new norm after orthonormalization'
1409 79 : call wrtout(std_out,message,'COLL')
1410 79 : call print_matlu(norm1%matlu(:),natom,prtopt=1)
1411 : end if
1412 :
1413 : ! ======================================================================
1414 : ! == Check that norm-identity is zero
1415 : ! ======================================================================
1416 79 : call init_oper(paw_dmft,norm2,opt_ksloc=2)
1417 79 : call init_oper(paw_dmft,norm3,opt_ksloc=2)
1418 79 : call identity_oper(norm2,2)
1419 79 : call add_matlu(norm1%matlu(:),norm2%matlu(:),norm3%matlu(:),natom,-1)
1420 79 : call destroy_oper(norm2)
1421 : if (pawprtvol > 2) then
1422 79 : write(message,'(2a)') ch10,' - Print new norm minus Identity '
1423 79 : call wrtout(std_out,message,'COLL')
1424 79 : call print_matlu(norm3%matlu(:),natom,prtopt=1,opt_exp=1)
1425 : end if
1426 79 : call destroy_oper(norm3)
1427 :
1428 79 : call destroy_oper(norm1)
1429 : ! call flush(std_out) ! debug debug debug debug
1430 : ! ABI_ERROR("Stop for debugging")
1431 :
1432 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1433 : ! New implementation, several atoms, general case.
1434 : else if (present(jkpt)) then !.or.iortho==2) then
1435 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1436 :
1437 91 : dimoverlap = 0
1438 182 : do iatom=1,natom
1439 91 : lpawu = paw_dmft%lpawu(iatom)
1440 91 : if (lpawu == -1) cycle
1441 182 : dimoverlap = dimoverlap + 2*lpawu + 1
1442 : ! write(6,*) "atom, dimoverlap",iatom,dimoverlap,natomcor
1443 : end do ! iatom
1444 :
1445 91 : dimoverlap = dimoverlap * nspinor
1446 :
1447 364 : ABI_MALLOC(largeoverlap,(dimoverlap,dimoverlap))
1448 364 : ABI_MALLOC(chipsivect,(dimoverlap,mbandc))
1449 273 : ABI_MALLOC(mat_tmp,(dimoverlap,mbandc))
1450 :
1451 : ! Big loop over isppol
1452 262 : do isppol=1,nsppol
1453 1230 : do ib=1,mbandc
1454 1059 : itot = 0
1455 2289 : do iatom=1,natom
1456 1059 : lpawu = paw_dmft%lpawu(iatom)
1457 1059 : if (lpawu == -1) cycle
1458 1059 : ndim = nspinor * (2*lpawu+1)
1459 6354 : chipsivect(itot+1:itot+ndim,ib) = paw_dmft%chipsi(1:ndim,ib,jkpt,isppol,iatom)
1460 1059 : itot = itot + ndim
1461 : !if(itot>dimoverlap) write(std_out,*) "itot>ndim",itot,ndim
1462 : ! write(6,*) "ib,iatom,im,ispinor",ib,iatom,im,ispinor,jkpt
1463 : end do ! iatom
1464 : end do ! ib
1465 :
1466 :
1467 : ! Calculation of overlap
1468 : call abi_xgemm("n","c",dimoverlap,dimoverlap,mbandc,cone,chipsivect(:,:),dimoverlap,&
1469 171 : & chipsivect(:,:),dimoverlap,czero,largeoverlap(:,:),dimoverlap)
1470 :
1471 : ! largeoverlap=czero
1472 : ! do ib=1,mbandc
1473 : ! do itot=1,dimoverlap
1474 : ! do itot1=1,dimoverlap
1475 : ! largeoverlap(itot,itot1)=largeoverlap(itot,itot1)+ &
1476 : !& psichivect(ib,itot)*conjg(psichivect(ib,itot1))
1477 : ! enddo ! itot1
1478 : ! enddo ! itot
1479 : ! enddo ! ib
1480 :
1481 : ! Math: orthogonalization of overlap
1482 171 : write(tag,'(i5)') jkpt
1483 171 : write(std_out,'(a,5x,2a)') ch10,"Overlap for ikpt = ",adjustl(tag)
1484 1026 : do itot=1,dimoverlap
1485 1026 : write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1486 : end do
1487 171 : call invsqrt_matrix(largeoverlap(:,:),dimoverlap,dum)
1488 171 : write(std_out,'(a,5x,2a)') ch10,"Inverse square root of overlap for ikpt = ",adjustl(tag)
1489 1026 : do itot=1,dimoverlap
1490 1026 : write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1491 : end do
1492 171 : write(std_out,'(a,5x,2a)') ch10,"Inverse square root of overlap for ikpt = ",adjustl(tag)
1493 1026 : do itot=1,dimoverlap
1494 1026 : write(std_out,'(100e9.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1495 : end do
1496 :
1497 : call abi_xgemm("n","n",dimoverlap,mbandc,dimoverlap,cone,largeoverlap(:,:),dimoverlap,&
1498 171 : & chipsivect(:,:),dimoverlap,czero,mat_tmp(:,:),dimoverlap)
1499 :
1500 : ! do ib=1,mbandc
1501 : ! wanall=czero
1502 : ! do itot=1,dimoverlap
1503 : ! do itot1=1,dimoverlap
1504 : ! wanall(itot)= wanall(itot)+psichivect(ib,itot1)*sqrtmatinv(itot,itot1)
1505 : ! enddo ! itot1
1506 : ! write(std_out,'(3i3,2x,i3,2x,2e15.5,2x,2e15.5)') jkpt,isppol,ib,itot,psichivect(ib,itot),wanall(itot)
1507 : ! enddo ! itot
1508 : ! iatomcor=0
1509 : ! do itot=1,dimoverlap
1510 : ! psichivect(ib,itot)=wanall(itot)
1511 : ! enddo
1512 : ! do iatom=1,natom
1513 : ! if(paw_dmft%lpawu(iatom).ne.-1) then
1514 : ! ndim=2*paw_dmft%lpawu(iatom)+1
1515 : ! iatomcor=iatomcor+1
1516 : ! do im=1,ndim
1517 : ! do ispinor=1,nspinor
1518 : ! paw_dmft%psichi(isppol,jkpt,ib,ispinor,iatom,im)=wanall(iatomcor,isppol,ispinor,im)
1519 : ! end do ! ispinor
1520 : ! end do ! im
1521 : ! endif
1522 : ! enddo ! iatom
1523 : ! enddo ! ib
1524 :
1525 :
1526 : ! Calculation of overlap (check)
1527 : call abi_xgemm("n","c",dimoverlap,dimoverlap,mbandc,cone,mat_tmp(:,:),dimoverlap,&
1528 171 : & mat_tmp(:,:),dimoverlap,czero,largeoverlap(:,:),dimoverlap)
1529 :
1530 : ! largeoverlap=czero
1531 : ! do ib=1,mbandc
1532 : ! do itot=1,dimoverlap
1533 : ! do itot1=1,dimoverlap
1534 : ! largeoverlap(itot,itot1)=largeoverlap(itot,itot1)+ &
1535 : !& psichivect(ib,itot)*conjg(psichivect(ib,itot1))
1536 : ! enddo ! itot1
1537 : ! enddo ! itot
1538 : ! enddo ! ib
1539 :
1540 171 : write(std_out,'(a,5x,2a)') ch10,"New overlap for ikpt = ",adjustl(tag)
1541 1026 : do itot=1,dimoverlap
1542 1026 : write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1543 : end do
1544 :
1545 : ! chipsivect -> chipsi
1546 1492 : do ib=1,mbandc
1547 1059 : itot = 0
1548 2289 : do iatom=1,natom
1549 1059 : lpawu = paw_dmft%lpawu(iatom)
1550 1059 : if (lpawu == -1) cycle
1551 1059 : ndim = nspinor * (2*lpawu+1)
1552 6354 : paw_dmft%chipsi(1:ndim,ib,jkpt,isppol,iatom) = mat_tmp(itot+1:itot+ndim,ib)
1553 1059 : itot = itot + ndim
1554 : end do ! iatom
1555 : end do ! ib
1556 :
1557 : ! End big loop over isppol
1558 : end do !isppol
1559 :
1560 91 : ABI_FREE(chipsivect)
1561 91 : ABI_FREE(largeoverlap)
1562 91 : ABI_FREE(mat_tmp)
1563 :
1564 : end if ! option
1565 :
1566 170 : end subroutine normalizechipsi
1567 : !!***
1568 :
1569 : !!****f* chipsi_renormalization/chipsi_gather
1570 : !! NAME
1571 : !! chipsi_gather
1572 : !!
1573 : !! FUNCTION
1574 : !! Gather chipsi from every CPU (parallelization over kpts).
1575 : !!
1576 : !! INPUTS
1577 : !! paw_dmft = data for DFT+DMFT calculations.
1578 : !!
1579 : !! OUTPUT
1580 : !!
1581 : !! NOTES
1582 : !!
1583 : !! SOURCE
1584 :
1585 96 : subroutine chipsi_gather(paw_dmft)
1586 :
1587 : !Arguments ------------------------------------
1588 : type(paw_dmft_type), intent(inout) :: paw_dmft
1589 : !Local variables ------------------------------
1590 : integer :: iatom,ib,ibuf,ierr,ikpt,irank,isppol,lpawu,mbandc,me_kpt
1591 : integer :: mkmem,natom,ndim,nproc,nspinor,nsppol,shift,siz_buf
1592 96 : integer, allocatable :: displs(:),recvcounts(:)
1593 96 : complex(dp), allocatable :: buffer(:),buffer_tot(:)
1594 : !************************************************************************
1595 :
1596 96 : me_kpt = paw_dmft%distrib%me_kpt
1597 96 : mkmem = paw_dmft%distrib%nkpt_mem(me_kpt+1)
1598 96 : mbandc = paw_dmft%mbandc
1599 96 : natom = paw_dmft%natom
1600 96 : nproc = paw_dmft%nproc
1601 96 : nspinor = paw_dmft%nspinor
1602 96 : nsppol = paw_dmft%nsppol
1603 96 : shift = paw_dmft%distrib%shiftk
1604 :
1605 96 : siz_buf = 0
1606 419 : do iatom=1,natom
1607 323 : lpawu = paw_dmft%lpawu(iatom)
1608 323 : if (lpawu == -1) cycle
1609 419 : siz_buf = siz_buf + 2*lpawu + 1
1610 : end do ! iatom
1611 96 : siz_buf = siz_buf * nspinor * mbandc * nsppol
1612 :
1613 288 : ABI_MALLOC(recvcounts,(nproc))
1614 192 : ABI_MALLOC(displs,(nproc))
1615 :
1616 408 : recvcounts(:) = siz_buf * paw_dmft%distrib%nkpt_mem(:)
1617 96 : displs(1) = 0
1618 312 : do irank=2,nproc
1619 312 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
1620 : end do ! irank
1621 :
1622 288 : ABI_MALLOC(buffer,(recvcounts(me_kpt+1)))
1623 288 : ABI_MALLOC(buffer_tot,(recvcounts(nproc)+displs(nproc)))
1624 :
1625 96 : ibuf = 0
1626 408 : do ikpt=1,mkmem
1627 1245 : do iatom=1,natom
1628 837 : lpawu = paw_dmft%lpawu(iatom)
1629 837 : if (lpawu == -1) cycle
1630 326 : ndim = nspinor * (2*lpawu+1)
1631 1147 : do isppol=1,nsppol
1632 6303 : do ib=1,mbandc
1633 30878 : buffer(ibuf+1:ibuf+ndim) = paw_dmft%chipsi(1:ndim,ib,ikpt+shift,isppol,iatom)
1634 5466 : ibuf = ibuf + ndim
1635 : end do ! ib
1636 : end do ! isppol
1637 : end do ! iatom
1638 : end do ! ikpt
1639 :
1640 96 : call xmpi_allgatherv(buffer(:),recvcounts(me_kpt+1),buffer_tot(:),recvcounts(:),displs(:),paw_dmft%distrib%comm_kpt,ierr)
1641 :
1642 96 : ibuf = 0
1643 849 : do ikpt=1,paw_dmft%nkpt
1644 2781 : do iatom=1,natom
1645 1932 : lpawu = paw_dmft%lpawu(iatom)
1646 1932 : if (lpawu == -1) cycle
1647 785 : ndim = nspinor * (2*lpawu+1)
1648 2722 : do isppol=1,nsppol
1649 14982 : do ib=1,mbandc
1650 70292 : paw_dmft%chipsi(1:ndim,ib,ikpt,isppol,iatom) = buffer_tot(ibuf+1:ibuf+ndim)
1651 13050 : ibuf = ibuf + ndim
1652 : end do ! ib
1653 : end do ! isppol
1654 : end do ! iatom
1655 : end do ! ikpt
1656 :
1657 96 : ABI_FREE(recvcounts)
1658 96 : ABI_FREE(displs)
1659 96 : ABI_FREE(buffer)
1660 96 : ABI_FREE(buffer_tot)
1661 :
1662 96 : end subroutine chipsi_gather
1663 : !!***
1664 :
1665 : end subroutine chipsi_renormalization
1666 : !!***
1667 :
1668 : !!****f* m_datafordmft/hybridization_asymptotic_coefficient
1669 : !! NAME
1670 : !! hybridization_asymptotic_coefficient
1671 : !!
1672 : !! FUNCTION
1673 : !! Compute some components for the limit of hybridization
1674 : !!
1675 : !! INPUTS
1676 : !! cryst_struc <type(crystal_t)>=crystal structure data
1677 : !! dft_occup
1678 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
1679 : !! pawtab <type(pawtab)>
1680 : !!
1681 : !! OUTPUT
1682 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
1683 : !!
1684 : !! NOTES
1685 : !!
1686 : !! SOURCE
1687 :
1688 34 : subroutine hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,hybri_coeff)
1689 :
1690 : !Arguments ------------------------------------
1691 : !scalars
1692 : type(crystal_t),intent(in) :: cryst_struc
1693 : type(paw_dmft_type), intent(in) :: paw_dmft
1694 : !type(pawang_type), intent(in) :: pawang
1695 : type(matlu_type), intent(inout) :: hybri_coeff(paw_dmft%natom)
1696 : !Local variables ------------------------------
1697 34 : type(oper_type) :: ham_a
1698 34 : type(oper_type) :: ham_b
1699 34 : type(oper_type) :: ham_squarelocal
1700 34 : type(oper_type) :: ham_squareks
1701 : integer :: iband1,iband2,ikpt,isppol
1702 : !************************************************************************
1703 :
1704 : ! call init_oper(paw_dmft,self_minus_hdc)
1705 34 : call init_oper(paw_dmft,ham_a)
1706 34 : call init_oper(paw_dmft,ham_b)
1707 34 : call init_oper(paw_dmft,ham_squareks)
1708 34 : call init_oper(paw_dmft,ham_squarelocal)
1709 :
1710 : ! Create self_minus_hdc%matlu = Sigma-Hdc in local orbitals
1711 : ! call add_matlu(self%oper(paw_dmft%dmft_nwlo)%matlu,self%hdc%matlu,&
1712 : !& self_minus_hdc%matlu,cryst_struc%natom,-1)
1713 :
1714 : ! write(message,'(a,2x,a)') ch10, " == self_minus_hdc (1)"
1715 : ! call wrtout(std_out,message,'COLL')
1716 : ! call print_matlu(self_minus_hdc%matlu,paw_dmft%natom,1,opt_exp=1)
1717 :
1718 : !! Create self_minus_hdc%ks = Upfold Sigma-Hdc
1719 : ! call upfold_oper(self_minus_hdc,paw_dmft,1)
1720 : ! call loc_oper(self_minus_hdc,paw_dmft,1)
1721 :
1722 : ! write(message,'(a,2x,a)') ch10, " == self_minus_hdc (2)"
1723 : ! call wrtout(std_out,message,'COLL')
1724 : ! call print_matlu(self_minus_hdc%matlu,paw_dmft%natom,1,opt_exp=1)
1725 :
1726 : ! Create ham_a%ks = H_ks in KS basis
1727 : !----------------------------------------------------
1728 333 : do iband1=1,paw_dmft%mbandc
1729 3752 : do iband2=1,paw_dmft%mbandc
1730 42588 : do ikpt=1,paw_dmft%nkpt
1731 91149 : do isppol=1,paw_dmft%nsppol
1732 87730 : if(iband1==iband2) then
1733 5540 : ham_a%ks(iband1,iband2,ikpt,isppol) = cmplx(paw_dmft%eigen_dft(iband1,ikpt,isppol),0.d0,kind=dp)
1734 : else
1735 43320 : ham_a%ks(iband1,iband2,ikpt,isppol) = czero
1736 : end if
1737 : end do
1738 : end do
1739 : end do
1740 : end do
1741 :
1742 : ! Create ham_a%matlu = H_ks in local orbitals
1743 : !---------------------------------------------
1744 34 : call downfold_oper(ham_a,paw_dmft)
1745 :
1746 : ! Symetrise the local quantity (energy levels)
1747 : !---------------------------------------------
1748 34 : call sym_matlu(ham_a%matlu,paw_dmft)
1749 :
1750 : ! Create ham_b%ks : Duplicate both ks and local part of ham_a into ham_b
1751 : !-----------------------------------------------------------------------
1752 34 : call copy_oper(ham_a,ham_b)
1753 :
1754 : ! Compute ham_squareks%ks : Make a product of the two KS version
1755 : !------------------------------------------------------------------
1756 34 : call prod_oper(ham_a,ham_b,ham_squareks,1)
1757 :
1758 : ! Compute ham_squareks%matlu
1759 : !---------------------------
1760 34 : call downfold_oper(ham_squareks,paw_dmft)
1761 :
1762 : ! Symetrise ham_squareks%matlu
1763 : !------------------------------
1764 34 : call sym_matlu(ham_squareks%matlu(:),paw_dmft)
1765 :
1766 : ! write(message,'(a,2x,a)') ch10, " == squareks"
1767 : ! call wrtout(std_out,message,'COLL')
1768 : ! call print_matlu(ham_squareks%matlu,paw_dmft%natom,1,opt_exp=1)
1769 :
1770 :
1771 : ! Compute ham_squarelocal%matlu
1772 : !-------------------------------
1773 34 : call prod_oper(ham_a,ham_b,ham_squarelocal,2)
1774 :
1775 : ! Compute the product in local orbitals
1776 : !--------------------------------------
1777 34 : call sym_matlu(ham_squarelocal%matlu(:),paw_dmft)
1778 :
1779 : ! write(message,'(a,2x,a)') ch10, " == squarelocal"
1780 : ! call wrtout(std_out,message,'COLL')
1781 : ! call print_matlu(ham_squarelocal%matlu,paw_dmft%natom,1,opt_exp=1)
1782 :
1783 : ! The difference of ham_squareks and ham_squarelocal
1784 : ! gives the coefficient that we are looking for ( such that F_ij(iw_n) = -C_ij/(iw_n) ).
1785 : !----------------------------------------------------------------------------------------
1786 34 : call add_matlu(ham_squareks%matlu(:),ham_squarelocal%matlu(:),hybri_coeff,cryst_struc%natom,-1)
1787 :
1788 : ! write(message,'(a,2x,a)') ch10, " == Coeff C_ij before sym"
1789 : ! call wrtout(std_out,message,'COLL')
1790 : ! call print_matlu(hybri_coeff,paw_dmft%natom,1,opt_exp=1)
1791 :
1792 : ! Symetrise the local quantity
1793 : !------------------------------
1794 34 : call sym_matlu(hybri_coeff,paw_dmft)
1795 :
1796 : ! write(message,'(a,2x,a)') ch10, " == Coeff C_ij after sym"
1797 : ! call wrtout(std_out,message,'COLL')
1798 : ! call print_matlu(hybri_coeff,paw_dmft%natom,1,opt_exp=1)
1799 :
1800 34 : call destroy_oper(ham_squarelocal)
1801 : ! call destroy_oper(self_minus_hdc)
1802 34 : call destroy_oper(ham_a)
1803 34 : call destroy_oper(ham_b)
1804 34 : call destroy_oper(ham_squareks)
1805 :
1806 :
1807 34 : end subroutine hybridization_asymptotic_coefficient
1808 : !!***
1809 :
1810 : !!****f* m_datafordmft/compute_wannier
1811 : !! NAME
1812 : !! compute_wannier
1813 : !!
1814 : !! FUNCTION
1815 : !! Compute the projected Wannier function in real space.
1816 : !!
1817 : !! INPUTS
1818 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
1819 : !! mpi_enreg=information about MPI parallelization
1820 : !!
1821 : !! OUTPUT
1822 : !!
1823 : !! NOTES
1824 : !!
1825 : !! SOURCE
1826 :
1827 0 : subroutine compute_wannier(paw_dmft,mpi_enreg)
1828 :
1829 : !Arguments ------------------------------------
1830 : type(paw_dmft_type), intent(inout) :: paw_dmft
1831 : type(MPI_type), intent(in) :: mpi_enreg
1832 : !Local variables ------------------------------
1833 : integer :: iatom,ibuf_psi,ibandc,ierr,iflavor,ik,ikpt,im,ispinor,isppol
1834 : integer :: itypat,lpawu,natom,nband_k,ndim,nkpt,nspinor,nsppol,siz_wan
1835 : !************************************************************************
1836 :
1837 0 : natom = paw_dmft%natom
1838 0 : nkpt = paw_dmft%nkpt
1839 0 : nspinor = paw_dmft%nspinor
1840 0 : nsppol = paw_dmft%nsppol
1841 :
1842 0 : ABI_MALLOC(paw_dmft%wannier,(paw_dmft%maxmeshsize,nspinor*(2*paw_dmft%maxlpawu+1)*nsppol,natom))
1843 0 : paw_dmft%wannier(:,:,:) = czero
1844 0 : ibuf_psi = 0
1845 :
1846 0 : do isppol=1,nsppol
1847 :
1848 0 : if (mpi_enreg%my_isppoltab(isppol) == 0) cycle
1849 0 : ik = 0
1850 :
1851 0 : do ikpt=1,nkpt
1852 :
1853 0 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
1854 0 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt)) cycle
1855 :
1856 0 : ik = ik + 1
1857 :
1858 0 : do ibandc=1,paw_dmft%mbandc
1859 0 : do ispinor=1,nspinor
1860 0 : do iatom=1,natom
1861 :
1862 0 : lpawu = paw_dmft%lpawu(iatom)
1863 0 : if (lpawu == -1) cycle
1864 0 : itypat = paw_dmft%typat(iatom)
1865 0 : ndim = 2*lpawu + 1
1866 0 : siz_wan = paw_dmft%radgrid(itypat)%mesh_size
1867 :
1868 0 : do im=1,ndim
1869 :
1870 0 : iflavor = im + (isppol*ispinor-1)*ndim
1871 :
1872 : paw_dmft%wannier(1:siz_wan,iflavor,iatom) = paw_dmft%wannier(1:siz_wan,iflavor,iatom) + &
1873 : & conjg(paw_dmft%chipsi(im+(ispinor-1)*ndim,ibandc,ikpt,isppol,iatom))* &
1874 0 : & paw_dmft%buf_psi(ibuf_psi+1:ibuf_psi+siz_wan)*paw_dmft%wtk(ikpt)
1875 :
1876 0 : ibuf_psi = ibuf_psi + siz_wan
1877 :
1878 : end do ! im
1879 :
1880 : end do ! iatom
1881 : end do ! ispinor
1882 : end do ! ibandc
1883 :
1884 : end do ! ikpt
1885 :
1886 : end do ! isppol
1887 :
1888 : ! No need to broadcast on every CPU
1889 0 : if (mpi_enreg%paral_kgb == 1 .and. mpi_enreg%nproc_band > 1) then
1890 0 : call xmpi_sum_master(paw_dmft%wannier(:,:,:),0,mpi_enreg%comm_band,ierr)
1891 : end if
1892 0 : call xmpi_sum_master(paw_dmft%wannier(:,:,:),0,mpi_enreg%comm_kpt,ierr)
1893 :
1894 0 : ABI_FREE(paw_dmft%buf_psi)
1895 :
1896 0 : end subroutine compute_wannier
1897 : !!***
1898 :
1899 : !!****f* m_datafordmft/print_wannier
1900 : !! NAME
1901 : !! print_wannier
1902 : !!
1903 : !! FUNCTION
1904 : !! Write projected Wannier functions on file.
1905 : !!
1906 : !! INPUTS
1907 : !! paw_dmft = data for self-consistent DFT+DMFT calculations.
1908 : !! istep = iteration step
1909 : !!
1910 : !! OUTPUT
1911 : !!
1912 : !! NOTES
1913 : !!
1914 : !! SOURCE
1915 :
1916 0 : subroutine print_wannier(paw_dmft,istep)
1917 :
1918 : !Arguments ------------------------------------
1919 : integer, intent(in) :: istep
1920 : type(paw_dmft_type), intent(in) :: paw_dmft
1921 : !Local variables ------------------------------
1922 : integer :: iatom,iflavor,ir,itypat,lpawu,nflavor,unt
1923 : character(len=3) :: tag_iter
1924 : character(len=10) :: tag_at
1925 : character(len=500) :: message
1926 : character(len=fnlen) :: tmpfil
1927 : !************************************************************************
1928 :
1929 0 : if (istep < 10) then
1930 0 : write(tag_iter,'("00",i1)') istep
1931 0 : else if (istep >= 10 .and. istep < 100) then
1932 0 : write(tag_iter,'("0",i2)') istep
1933 0 : else if (istep >= 100 .and. istep < 1000) then
1934 0 : write(tag_iter,'(i3)') istep
1935 : else
1936 0 : tag_iter="xxx"
1937 : end if ! istep
1938 :
1939 0 : do iatom=1,paw_dmft%natom
1940 0 : lpawu = paw_dmft%lpawu(iatom)
1941 0 : if (lpawu == -1) cycle
1942 0 : itypat = paw_dmft%typat(iatom)
1943 0 : nflavor = (2*lpawu+1) * paw_dmft%nspinor * paw_dmft%nsppol
1944 0 : call int2char4(iatom,tag_at)
1945 0 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
1946 0 : tmpfil = trim(paw_dmft%filapp)//'_Wannier_functions_iatom'//trim(tag_at)//'_'//tag_iter
1947 0 : if (open_file(tmpfil,message,newunit=unt) /= 0) ABI_ERROR(message)
1948 0 : write(unt,'(4a)') "# Radial part of projective Wannier functions, after orthonormalization, for each flavor.", &
1949 0 : & ch10,"# First column is the radius (bohr) and the other columns correspond to the radial part of flavor i,", &
1950 0 : & " where i=1...2*(2*l+1) (spins up are first)"
1951 0 : do ir=1,paw_dmft%radgrid(itypat)%mesh_size
1952 0 : write(unt,*) paw_dmft%radgrid(itypat)%rad(ir),(dble(paw_dmft%wannier(ir,iflavor,iatom)),iflavor=1,nflavor)
1953 : end do ! ir
1954 0 : close(unt)
1955 : end do ! iatom
1956 :
1957 0 : end subroutine print_wannier
1958 : !!***
1959 :
1960 : END MODULE m_datafordmft
1961 : !!***
|
