Branch data 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 : 101 : subroutine datafordmft(cg,cprj,cryst_struc,dft_occup,dimcprj,dtset,eigen,mband_cprj,mcg,&
105 : 101 : & 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 : 101 : type(oper_type) :: loc_norm_check
134 : 101 : integer, allocatable :: displs(:),recvcounts(:)
135 : 101 : complex(dp), allocatable :: buf_chipsi(:),buf_chipsi_tot(:),chipsi_tmp(:),cwprj(:,:)
136 : 101 : type(pawcprj_type), allocatable :: cwaveprj(:,:)
137 : 101 : 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 : 101 : mband = paw_dmft%mband
144 : 101 : mbandc = paw_dmft%mbandc
145 : 101 : mkmem = paw_dmft%mkmem
146 : 101 : natom = paw_dmft%natom
147 : 101 : nkpt = paw_dmft%nkpt
148 : 101 : nspinor = paw_dmft%nspinor
149 : 101 : nsppol = paw_dmft%nsppol
150 [ + + ]: 240 : nsppol_mem = sum(mpi_enreg%my_isppoltab(1:nsppol))
151 : 101 : pawprtvol = dtset%pawprtvol
152 : 101 : prt_wan = (paw_dmft%dmft_prtwan == 1)
153 : :
154 [ + + ]: 101 : if (abs(pawprtvol) >= 3) then
155 : 77 : write(message,*) " number of k-points used is nkpt=nkpt ",nkpt
156 : 77 : call wrtout(std_out,message,'COLL')
157 : 77 : write(message,*) " warning: parallelized version ",nkpt
158 : 77 : call wrtout(std_out,message,'COLL')
159 : 77 : write(message,*) " weights k-points used is wtk=wtk"
160 : 77 : call wrtout(std_out,message,'COLL')
161 : : end if ! abs(pawprtvol)>=3
162 : :
163 [ - + ]: 101 : if (usecprj == 0) then
164 : 0 : write(message,*) " usecprj=0 : BUG in datafordmft",usecprj
165 : 0 : ABI_BUG(message)
166 : : end if
167 : :
168 [ - + ]: 101 : 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 : 101 : paral_kgb = mpi_enreg%paral_kgb
182 [ + + ]: 101 : comm_kpt = merge(mpi_enreg%comm_kpt,mpi_enreg%comm_cell,paral_kgb==1)
183 : 101 : comm_band = mpi_enreg%comm_band
184 : 101 : me_kpt = mpi_enreg%me_kpt
185 : 101 : me_band = mpi_enreg%me_band
186 : 101 : nproc_band = mpi_enreg%nproc_band
187 : 101 : nproc_spkpt = mpi_enreg%nproc_spkpt
188 : :
189 [ - + ]: 101 : 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 : 101 : iorder_cprj = 0
195 [ - + ]: 101 : 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 : 101 : nbandi = paw_dmft%dmftbandi
202 : 101 : nbandf = paw_dmft%dmftbandf
203 : 101 : t2g = (paw_dmft%dmft_t2g == 1)
204 : 101 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
205 : : !use_full_chipsi = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
206 [ + - ]: 101 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
207 : : use_full_chipsi = .true.
208 : : else
209 : 101 : use_full_chipsi = (paw_dmft%dmft_full_chipsi /= 0)
210 : : endif
211 : :
212 [ - + - - ]: 101 : 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 : 101 : write(message,'(2a)') ch10,' == Prepare data for DFT+DMFT calculation'
221 : 101 : call wrtout(std_out,message,'COLL')
222 [ + + ]: 101 : if (abs(pawprtvol) >= 3) then
223 : 77 : write(message,'(2a)') ch10,'---------------------------------------------------------------'
224 : : ! call wrtout(ab_out,message,'COLL')
225 : 77 : call wrtout(std_out,message,'COLL')
226 : 77 : write(message,'(12a)') ch10,' Print useful data (as a check)',ch10,&
227 : 77 : & ' - Overlap of KS wfc with atomic orbital inside sphere',ch10,&
228 : 77 : & ' - Eigenvalues',ch10,&
229 : 77 : & ' - Weights of k-points',ch10,&
230 : 77 : & ' - Number of spins ',ch10,&
231 : 154 : & ' - Number of states'
232 : : ! call wrtout(ab_out,message,'COLL')
233 : 77 : call wrtout(std_out,message,'COLL')
234 : 77 : write(message,'(2a)') ch10,'---------------------------------------------------------------'
235 : 77 : call wrtout(std_out,message,'COLL')
236 : : end if ! abs(pawprtvol)>=3
237 : :
238 [ - + - - ]: 101 : 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 : 101 : maxlpawu = paw_dmft%maxlpawu
245 : 101 : maxmeshsize = paw_dmft%maxmeshsize
246 : 101 : maxnproju = paw_dmft%maxnproju
247 : :
248 : : !***************** in forlb.eig
249 [ + + + + ]: 101 : if (paw_dmft%myproc == 0 .and. abs(pawprtvol) >= 3) then
250 [ - + ]: 26 : if (open_file('forlb.eig',message,newunit=unt,form='formatted',status='unknown') /= 0) ABI_ERROR(message)
251 : 26 : rewind(unt)
252 : 26 : write(unt,*) "Number of bands, spins, and k-point; and spin-orbit flag"
253 : 26 : write(unt,*) mband,nsppol,nkpt,my_nspinor,nbandi,nbandf
254 : 26 : write(unt,*) " For each k-point, eigenvalues for each band"
255 [ + - + + ]: 275 : write(unt,*) (dtset%wtk(ikpt),ikpt=1,nkpt)
256 : 26 : band_index = 0
257 [ + + ]: 57 : do isppol=1,nsppol
258 : 31 : write(unt,*) " For spin"
259 : 31 : write(unt,*) isppol
260 [ + + ]: 396 : do ikpt=1,nkpt
261 : 339 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
262 : 339 : ibandc = 0
263 : 339 : write(unt,*) " For k-point"
264 : 339 : write(unt,*) ikpt
265 [ + + ]: 9427 : do iband=1,mband
266 [ + + ]: 9427 : if (paw_dmft%band_in(iband)) then
267 : 4363 : ibandc = ibandc + 1
268 : 4363 : write(unt,'(2i6,4x,f20.15)') ibandc,ikpt,eigen(iband+band_index)*two
269 : : end if
270 : : end do ! iband
271 : 370 : band_index = band_index + nband_k
272 : : end do ! ikpt
273 : : end do ! isppol
274 : 26 : 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 [ + + ]: 101 : 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 : 101 : paw_dmft%eigen => eigen(:)
294 : 101 : band_index = 0
295 [ + + ]: 240 : do isppol=1,nsppol
296 [ + + ]: 1426 : do ikpt=1,nkpt
297 : 1186 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
298 : 1186 : ibandc = 0
299 [ + + ]: 35042 : do iband=1,mband
300 [ + + ]: 35042 : if (paw_dmft%band_in(iband)) then
301 : 15688 : 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 [ + + ]: 15688 : if (isppol .eq. 1) then
305 : 12994 : 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 : 1325 : band_index = band_index + nband_k
312 : : end do ! ikpt
313 : : end do ! isppol
314 : :
315 [ + + ]: 101 : if (abs(pawprtvol) >= 3) then
316 : 77 : write(message,'(2a)') ch10,' datafordmft : eigenvalues written on file'
317 : 77 : 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 [ - + + - : 868 : 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 : 101 : call pawcprj_alloc(cwaveprj(:,:),0,dimcprj(:))
330 : : !write(std_out,*) "after alloc cprj"
331 : :
332 [ - + + - : 404 : ABI_MALLOC(cwprj,(maxnproju,2*maxlpawu+1))
+ - - + +
- + - + -
- + - + ]
333 [ - + + - : 303 : ABI_MALLOC(chipsi_tmp,(maxmeshsize))
+ - + - -
+ - + ]
334 : :
335 : 101 : siz_buf = 0
336 : 101 : siz_buf_psi = 0
337 [ + + ]: 442 : do iatom=1,natom
338 : 341 : lpawu = paw_dmft%lpawu(iatom)
339 [ + + ]: 341 : if (lpawu == -1) cycle
340 : 121 : itypat = paw_dmft%typat(iatom)
341 : 121 : ndim = 2*lpawu + 1
342 : 121 : siz_buf = siz_buf + ndim
343 : 442 : siz_buf_psi = siz_buf_psi + ndim*paw_dmft%radgrid(itypat)%mesh_size
344 : : end do ! iatom
345 : 101 : siz_buf = siz_buf * mbandc * nspinor * mkmem * nsppol_mem
346 : 101 : siz_buf_psi = siz_buf_psi * mbandc * nspinor * mkmem * nsppol_mem
347 : :
348 [ - + + - : 303 : ABI_MALLOC(recvcounts,(nproc_spkpt))
+ - + - -
+ - + ]
349 [ - + + - : 202 : ABI_MALLOC(displs,(nproc_spkpt))
- + - + ]
350 : 101 : call xmpi_allgather(siz_buf,recvcounts(:),comm_kpt,ierr)
351 : :
352 : 101 : displs(1) = 0
353 [ + + ]: 209 : do irank=2,nproc_spkpt
354 : 209 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
355 : : end do ! irank
356 : :
357 [ - + + - : 303 : ABI_MALLOC(buf_chipsi,(siz_buf))
+ - + - -
+ - + ]
358 [ - + + - : 303 : ABI_MALLOC(buf_chipsi_tot,(recvcounts(nproc_spkpt)+displs(nproc_spkpt)))
+ - + - -
+ - + ]
359 [ - + ]: 101 : if (prt_wan) then
360 [ # # # # : 0 : ABI_MALLOC(paw_dmft%buf_psi,(siz_buf_psi))
# # # # #
# # # #
# ]
361 : : end if
362 : :
363 : 101 : icprj = 0
364 : 101 : icg = 0
365 : 101 : ibuf_psi = 0
366 : 101 : ibuf_chipsi = 0
367 : :
368 [ + + ]: 88692 : buf_chipsi(:) = czero
369 : :
370 [ + + ]: 240 : do isppol=1,nsppol
371 : :
372 [ + + ]: 139 : if (mpi_enreg%my_isppoltab(isppol) == 0) cycle
373 : 119 : ik = 0
374 : :
375 [ + + ]: 1142 : do ikpt=1,nkpt
376 : :
377 : 922 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
378 [ + + ]: 922 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle
379 : :
380 : 649 : nband_k_cprj = nband_k / nproc_band
381 : 649 : npw = paw_dmft%npwarr(ikpt)
382 : 649 : icgb = icg
383 : 649 : ik = ik + 1
384 : 649 : ibandc = 0
385 : : ! LOOP OVER BANDS
386 : 649 : ib = 0
387 : :
388 [ + + ]: 18449 : do iband=1,nband_k
389 : :
390 : : ! Parallelization: treat only some bands
391 : 17800 : verif_band = .true.
392 [ + + ]: 17800 : if (paral_kgb == 1) then
393 [ + + ]: 6400 : if (mod((iband-1)/mpi_enreg%bandpp,nproc_band) /= me_band) verif_band = .false.
394 : : else
395 [ + - ]: 11400 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol) /= me_kpt) verif_band = .false.
396 : : end if
397 : :
398 : 13000 : if (verif_band) ib = ib + 1
399 : :
400 [ + + ]: 17800 : if (paw_dmft%band_in(iband)) then
401 : 7435 : ibandc = ibandc + 1
402 : : else
403 : 10365 : icgb = icgb + npw*nspinor
404 : 10365 : cycle
405 : : end if
406 : :
407 [ + + ]: 7435 : 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 : 5635 : & proc_distrb=mpi_enreg%proc_distrb(:,:,:))
412 : : end if
413 : :
414 [ + + ]: 16943 : do ispinor=1,my_nspinor
415 [ + + ]: 39414 : do iatom=1,natom
416 : 30555 : lpawu = paw_dmft%lpawu(iatom)
417 [ + + ]: 30555 : if (lpawu == -1) cycle
418 [ + + ]: 15011 : lpawu1 = merge(2,lpawu,t2g.or.x2my2d)
419 : 15011 : itypat = paw_dmft%typat(iatom)
420 : 15011 : ndim = 2*lpawu + 1
421 : 15011 : nproju = pawtab(itypat)%nproju
422 : 15011 : siz_proj = paw_dmft%siz_proj(itypat)
423 : 15011 : siz_wan = paw_dmft%radgrid(itypat)%mesh_size
424 : 15011 : siz_paw = min(siz_wan,paw_dmft%int_meshsz(itypat))
425 : 15011 : rint = paw_dmft%radgrid(itypat)%rad(siz_proj)
426 : :
427 [ + + ]: 15011 : if (verif_band) then
428 : 13211 : lmn_size = pawtab(itypat)%lmn_size
429 [ + + ]: 381601 : do ilmn=1,lmn_size
430 : : ! ------------ Select l=lpawu.
431 [ + + ]: 368390 : if (pawtab(itypat)%indlmn(1,ilmn) /= lpawu1) cycle
432 : :
433 : 162742 : im = pawtab(itypat)%indlmn(2,ilmn) + lpawu1 + 1
434 [ - + ]: 162742 : if (x2my2d) then
435 [ # # ]: 0 : if (im /= 5) cycle
436 : : im = 1
437 [ + + ]: 162742 : else if (t2g) then
438 [ + + ]: 10200 : if (im == 3 .or. im == 5) cycle
439 [ + + ]: 6120 : if (im == 4) im = 3
440 : : end if
441 : 158662 : iproj = pawtab(itypat)%indlmn(3,ilmn)
442 : : cwprj(iproj,im) = cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn), &
443 : 381601 : & cwaveprj(iatom,ispinor)%cp(2,ilmn),kind=dp)
444 : : end do ! ilmn
445 : : end if ! verif
446 : :
447 [ + + ]: 112461 : do im=1,ndim
448 : :
449 : 88591 : ibuf_chipsi = ibuf_chipsi + 1
450 : :
451 [ - + ]: 88591 : 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 [ + + + + ]: 247253 : 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 : 119146 : ibuf_psi = ibuf_psi + siz_wan
491 : :
492 : : end do ! im
493 : : end do ! iatom
494 : 16294 : icgb = icgb + npw
495 : : end do ! ispinor
496 : : end do ! iband
497 : :
498 : 649 : icprj = icprj + nband_k_cprj*nspinor
499 : 1061 : 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 : 101 : call pawcprj_free(cwaveprj)
521 [ + - + + : 459 : ABI_FREE(cwaveprj)
- + - + -
+ ]
522 : :
523 [ - + ]: 101 : ABI_FREE(cwprj)
524 [ - + ]: 101 : ABI_FREE(chipsi_tmp)
525 : :
526 : : !==========================================================================
527 : : !********************* Gather information for MPI before printing
528 : : !==========================================================================
529 : :
530 : : !call xmpi_barrier(comm_kpt)
531 [ + + ]: 101 : 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 : 101 : call xmpi_allgatherv(buf_chipsi(:),siz_buf,buf_chipsi_tot(:),recvcounts(:),displs(:),comm_kpt,ierr)
536 : :
537 : 101 : ABI_FREE(displs)
538 : 101 : ABI_FREE(recvcounts)
539 [ - + ]: 101 : ABI_FREE(buf_chipsi)
540 : :
541 : : ! Reorder the chipsi since the kpts can have any arbitrary distribution over the different MPI processes
542 : 101 : ibuf_chipsi = 0
543 [ + + ]: 310 : do irank=0,nproc_spkpt-1
544 [ + + ]: 617 : do isppol=1,nsppol
545 [ + + ]: 3850 : do ikpt=1,nkpt
546 : :
547 : 3334 : nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
548 [ + + ]: 3334 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,irank)) cycle
549 : :
550 [ + + ]: 17181 : do ibandc=1,mbandc
551 [ + + ]: 40214 : do ispinor=1,nspinor
552 [ + + ]: 107056 : do iatom=1,natom
553 : 70176 : lpawu = paw_dmft%lpawu(iatom)
554 [ + + ]: 70176 : if (lpawu == -1) cycle
555 : 44600 : ndim = 2*lpawu + 1
556 [ + + ]: 331656 : paw_dmft%chipsi(1+(ispinor-1)*ndim:ispinor*ndim,ibandc,ikpt,isppol,iatom) = buf_chipsi_tot(ibuf_chipsi+1:ibuf_chipsi+ndim)
557 : 46768 : 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 [ - + ]: 101 : 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 [ + + ]: 101 : if (paw_dmft%myproc == 0) then
591 : 44 : 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 [ + + ]: 101 : if (paw_dmft%dmft_kspectralfunc == 0) then
598 : :
599 : 97 : call init_oper(paw_dmft,loc_norm_check,opt_ksloc=2)
600 : 97 : 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 : 97 : write(message,'(2a)') ch10," == The DMFT orbitals are now projected on the correlated bands"
610 : 97 : call wrtout(std_out,message,'COLL')
611 : :
612 : 97 : write(message,'(2a,i4)') ch10," == Check: Downfolded Occupations and Norm of unnormalized projected orbitals"
613 : 97 : call wrtout(std_out,message,'COLL')
614 : :
615 [ + + ]: 97 : 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 : 97 : call sym_matlu(dft_occup%matlu(:),paw_dmft)
631 : :
632 : 97 : write(message,'(2a,i4)') ch10,' ------ Symmetrized Occupations'
633 : 97 : call wrtout(std_out,message,'COLL')
634 : :
635 : 97 : call print_matlu(dft_occup%matlu(:),natom,pawprtvol)
636 : :
637 : : ! symmetrize and print norms
638 : 97 : call sym_matlu(loc_norm_check%matlu(:),paw_dmft)
639 : :
640 : 97 : write(message,'(2a,i4)') ch10,' ------ Symmetrized Norm'
641 : 97 : call wrtout(std_out,message,'COLL')
642 : :
643 : 97 : call print_matlu(loc_norm_check%matlu(:),natom,pawprtvol)
644 : :
645 : : ! Tests density matrix DFT+U and density matrix computed here.
646 [ + - + + ]: 97 : 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 : 91 : write(message,'(2a)') ch10,&
694 : 182 : & ' Warning: Consistency of density matrices computed from projection has not been checked: use dmftcheck>=2 '
695 : 91 : call wrtout(std_out,message,'COLL')
696 : : end if
697 : :
698 : 97 : call destroy_oper(loc_norm_check)
699 : : end if ! dmft_kspectralfunc=0
700 : :
701 [ + - - + : 202 : if (present(nbandkss)) then
- - - - -
+ ]
702 [ + + + + : 101 : 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 : 97 : 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 : 97 : real(dp), allocatable :: occ_dft(:,:,:)
741 : : ! *********************************************************************
742 : :
743 [ - + + - : 485 : ABI_MALLOC(occ_dft,(mbandc,nkpt,nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
744 : :
745 : 97 : band_index = 0
746 [ + + ]: 232 : do isppol=1,nsppol
747 [ + + ]: 1374 : do ikpt=1,nkpt
748 : 1142 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
749 : 1142 : ibandc = 0
750 [ + + ]: 34118 : do iband=1,nband_k
751 [ + + ]: 34118 : if (paw_dmft%band_in(iband)) then
752 : 15292 : ibandc = ibandc + 1
753 : 15292 : occ_dft(ibandc,ikpt,isppol) = occ(iband+band_index)
754 : : end if
755 : : end do ! iband
756 : 1277 : band_index = band_index + nband_k
757 : : end do ! ikpt
758 : : end do ! isppol
759 : :
760 [ + + + + : 4945 : if (nsppol == 1 .and. my_nspinor == 1) occ_dft(:,:,:) = occ_dft(:,:,:) * half
+ + + + +
+ ]
761 : :
762 : 97 : call downfold_oper(xnorm_check,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=2)
763 : 97 : 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 : 97 : & iproc=paw_dmft%distrib%me_kpt,option=3,op_ks_diag=occ_dft(:,:,:))
766 : 97 : call xmpi_matlu(xocc_check%matlu(:),natom,paw_dmft%distrib%comm_kpt)
767 : :
768 [ - + ]: 97 : ABI_FREE(occ_dft)
769 : :
770 : 97 : 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 : 45 : 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 : 45 : t2g = (paw_dmft%dmft_t2g == 1)
806 : 45 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
807 : :
808 [ - + ]: 45 : if (open_file('forlb.ovlp',msg,newunit=unt,form='formatted',status='unknown') /= 0) ABI_ERROR(msg)
809 : 45 : rewind(unt)
810 : :
811 : : ! Header for calc_uCRPA.F90
812 : 45 : write(unt,*) "# isppol nspinor natom m Re(<chi|psi>) Im(<chi|psi>)"
813 [ + + + + : 135 : if (count(pawtab(:)%lpawu /= -1) == 1) then
+ - ]
814 [ + + ]: 135 : do itypat=1,paw_dmft%ntypat
815 : 90 : lpawu = pawtab(itypat)%lpawu
816 [ + + ]: 90 : if (lpawu == -1) cycle
817 [ + + ]: 90 : if (t2g) then
818 : 17 : write(unt,*) "l= ",1,itypat
819 [ - + ]: 28 : else if (x2my2d) then
820 : 0 : write(unt,*) "l= ",0,itypat
821 : : else
822 : 28 : 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 : 45 : write(unt,*) "Bands ",paw_dmft%dmftbandi,paw_dmft%dmftbandf
830 : :
831 [ + + ]: 113 : do isppol=1,paw_dmft%nsppol
832 [ + + ]: 623 : do ikpt=1,paw_dmft%nkpt
833 : : ! rewind(1023)
834 : 510 : write(unt,'(a6,2x,i6)') "ikpt =",ikpt
835 : 510 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*paw_dmft%nkpt)
836 : 510 : ibandc = 0
837 [ + + ]: 13798 : do iband=1,nband_k
838 [ + + ]: 13220 : if (paw_dmft%band_in(iband)) then
839 : 5734 : ibandc = ibandc + 1
840 : 5734 : write(unt,'(a8,2x,i6)') " iband =",iband
841 : : else
842 : : cycle
843 : : end if
844 [ + + ]: 13402 : do ispinor=1,paw_dmft%nspinor
845 [ + + ]: 42032 : do iatom=1,paw_dmft%natom
846 : 21654 : lpawu = paw_dmft%lpawu(iatom)
847 [ + + ]: 21654 : if (lpawu == -1) cycle
848 : 13310 : ndim = 2*lpawu + 1
849 [ + + ]: 104154 : do im=1,ndim
850 [ + + ]: 83686 : if (t2g) then
851 : 3060 : im1 = mt2g(im)
852 [ - + ]: 80626 : else if (x2my2d) then
853 : 0 : im1 = 5
854 : : else
855 : 80626 : im1 = im
856 : : end if
857 : :
858 : 83686 : write(unt,'(4i6,3x,2f23.15)') isppol,ispinor,iatom,im1,&
859 : 83686 : & dble(paw_dmft%chipsi(im+(ispinor-1)*ndim,ibandc,ikpt,isppol,iatom)),&
860 : 189026 : & 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 : 45 : close(unt)
870 : :
871 : 45 : 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 : 217 : 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 : 217 : type(matlu_type), allocatable :: levels_temp(:),magnfield(:)
906 : : !************************************************************************
907 : :
908 : 217 : natom = paw_dmft%natom
909 : 217 : nspinor = paw_dmft%nspinor
910 : 217 : nsppol = paw_dmft%nsppol
911 [ + + ]: 217 : 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 : 217 : & iproc=paw_dmft%distrib%me_kpt,op_ks_diag=paw_dmft%eigen_dft(:,:,:))
919 : 217 : 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 [ + + ]: 983 : do iatom=1,natom
924 : 766 : lpawu = paw_dmft%lpawu(iatom)
925 [ + + ]: 766 : if (lpawu == -1) cycle
926 : 281 : ndim = paw_dmft%nspinor * (2*lpawu+1)
927 [ + + ]: 707 : do isppol=1,paw_dmft%nsppol
928 [ + + ]: 2945 : do im=1,ndim
929 : 2664 : 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 [ + + + + : 18221 : 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 : 217 : call sym_matlu(energy_level%matlu(:),paw_dmft)
937 [ + + ]: 217 : if (present(nondiag)) then
938 : 148 : call checkdiag_matlu(energy_level%matlu(:),natom,tol7,nondiag)
939 : : end if
940 : :
941 : 217 : write(tag,'(f13.5)') paw_dmft%fermie
942 : 217 : write(message,'(a,2x,2a)') ch10," == Print Energy levels in cubic basis for Fermi Level = ",adjustl(tag)
943 : 217 : call wrtout(std_out,message,'COLL')
944 : : !call print_oper(energy_level,1,paw_dmft,1)
945 : 217 : call print_matlu(energy_level%matlu(:),natom,1)
946 : :
947 : : !========================================================
948 : : !Add Zeeman contibutions to local energy levels in slm
949 : : !========================================================
950 [ + + + - ]: 217 : 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 : 217 : 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 : 100 : 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 : 100 : type(oper_type) :: norm,oper_temp
1007 : : character(len=500) :: message
1008 : 100 : real(dp), allocatable :: wtk_tmp(:)
1009 : : character(len=5) :: tag
1010 : : ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:)
1011 : : !************************************************************************
1012 : :
1013 : : DBG_ENTER("COLL")
1014 : :
1015 : 100 : option = 3
1016 [ - + ]: 100 : if (present(opt)) then
1017 [ - + ]: 100 : if (opt == 2 .or. opt == 3) option = opt
1018 : : end if
1019 : 100 : pawprtvol = 2
1020 : :
1021 : 100 : natom = paw_dmft%natom
1022 : 100 : 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 : 100 : write(message,'(6a)') ch10, &
1053 : 100 : & ' =================================================== ',&
1054 : 100 : & ch10,' == The DMFT orbitals will now be orthonormalized == ',&
1055 : 200 : & ch10,' =================================================== '
1056 : 100 : call wrtout(std_out,message,'COLL')
1057 : :
1058 [ + + ]: 100 : if (option == 2) then ! option==2
1059 : : ! ====================================
1060 : : ! == renormalize k-point after k-point
1061 : : ! ====================================
1062 : :
1063 : 21 : write(tag,'(i5)') nkpt
1064 : 21 : write(message,'(3a)') ch10,' Number of k-points: ',adjustl(tag)
1065 : 21 : call wrtout(std_out,message,'COLL')
1066 [ + + ]: 389 : do jkpt=1,nkpt ! jkpt
1067 : 368 : write(tag,'(i5)') jkpt
1068 : 368 : write(message,'(3a)') ch10,' == Orthonormalization for k-point: ',adjustl(tag)
1069 : 368 : call wrtout(std_out,message,'COLL')
1070 [ + + ]: 368 : if (paw_dmft%distrib%procb(jkpt) /= paw_dmft%distrib%me_kpt) cycle
1071 : 389 : call normalizechipsi(1,paw_dmft,jkpt=jkpt)
1072 : : end do ! jkpt
1073 : 21 : write(message,'(2a)') ch10,' ===== Every k-point has been orthonormalized'
1074 : 21 : 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 : 100 : call chipsi_gather(paw_dmft)
1091 : :
1092 : : !== Change back repr for norm
1093 : :
1094 : : !===============================================
1095 : : !== Compute norm with new chipsi
1096 : : !===============================================
1097 : :
1098 : 100 : write(message,'(2a)') ch10,' ===== Compute new norm after renormalization'
1099 : 100 : call wrtout(std_out,message,'COLL')
1100 : 100 : call init_oper(paw_dmft,oper_temp,opt_ksloc=2)
1101 : 100 : call identity_oper(oper_temp,2)
1102 : :
1103 [ + + - + : 100 : 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 : 96 : write(message,'(2a)') ch10,' ===== Starting downfold'
1123 : 96 : call wrtout(std_out,message,'COLL')
1124 : 96 : call init_oper(paw_dmft,norm,opt_ksloc=2)
1125 : 96 : call downfold_oper(norm,paw_dmft,procb=paw_dmft%distrib%procb(:),iproc=paw_dmft%distrib%me_kpt,option=2)
1126 : 96 : call xmpi_matlu(norm%matlu(:),natom,paw_dmft%distrib%comm_kpt)
1127 : 96 : write(message,'(2a)') ch10,' ===== Finished downfold'
1128 : 96 : 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 : 96 : 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 : 96 : call diff_oper('Overlap after renormalization','Identity',norm,oper_temp,0,tol6)
1151 : :
1152 : : end if ! dmft_kspectralfunc
1153 : :
1154 : 100 : call destroy_oper(norm)
1155 : 100 : call destroy_oper(oper_temp)
1156 : :
1157 [ - + - - : 100 : 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 : 186 : 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 : 186 : type(oper_type) :: norm1,norm2,norm3
1192 : 186 : real(dp), allocatable :: si(:),sr(:)
1193 : 186 : 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 : 186 : mbandc = paw_dmft%mbandc
1205 : 186 : natom = paw_dmft%natom
1206 : 186 : nspinor = paw_dmft%nspinor
1207 : 186 : nsppol = paw_dmft%nsppol
1208 : 186 : ndim_max = (2*paw_dmft%maxlpawu+1) * nspinor
1209 : 186 : pawprtvol = 3
1210 : :
1211 [ + + - + ]: 186 : 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 [ + + ]: 114 : 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 : 107 : dimoverlap = 0
1438 [ + + ]: 262 : do iatom=1,natom
1439 : 155 : lpawu = paw_dmft%lpawu(iatom)
1440 [ + + ]: 155 : if (lpawu == -1) cycle
1441 : 262 : dimoverlap = dimoverlap + 2*lpawu + 1
1442 : : ! write(6,*) "atom, dimoverlap",iatom,dimoverlap,natomcor
1443 : : end do ! iatom
1444 : :
1445 : 107 : dimoverlap = dimoverlap * nspinor
1446 : :
1447 [ - + + - : 428 : ABI_MALLOC(largeoverlap,(dimoverlap,dimoverlap))
- + + - +
- + - - +
- + ]
1448 [ - + + - : 428 : ABI_MALLOC(chipsivect,(dimoverlap,mbandc))
- + + - +
- + - - +
- + ]
1449 [ - + + - : 321 : ABI_MALLOC(mat_tmp,(dimoverlap,mbandc))
- + + - -
+ - + ]
1450 : :
1451 : : ! Big loop over isppol
1452 [ + + ]: 294 : do isppol=1,nsppol
1453 [ + + ]: 2606 : do ib=1,mbandc
1454 : 2419 : itot = 0
1455 [ + + ]: 9105 : do iatom=1,natom
1456 : 6499 : lpawu = paw_dmft%lpawu(iatom)
1457 [ + + ]: 6499 : if (lpawu == -1) cycle
1458 : 5139 : ndim = nspinor * (2*lpawu+1)
1459 [ + + ]: 67554 : chipsivect(itot+1:itot+ndim,ib) = paw_dmft%chipsi(1:ndim,ib,jkpt,isppol,iatom)
1460 : 3779 : 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 : 187 : & 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 : 187 : write(tag,'(i5)') jkpt
1483 : 187 : write(std_out,'(a,5x,2a)') ch10,"Overlap for ikpt = ",adjustl(tag)
1484 [ + + ]: 1714 : do itot=1,dimoverlap
1485 : 1714 : write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1486 : : end do
1487 : 187 : call invsqrt_matrix(largeoverlap(:,:),dimoverlap,dum)
1488 : 187 : write(std_out,'(a,5x,2a)') ch10,"Inverse square root of overlap for ikpt = ",adjustl(tag)
1489 [ + + ]: 1714 : do itot=1,dimoverlap
1490 : 1714 : write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1491 : : end do
1492 : 187 : write(std_out,'(a,5x,2a)') ch10,"Inverse square root of overlap for ikpt = ",adjustl(tag)
1493 [ + + ]: 1714 : do itot=1,dimoverlap
1494 : 1714 : 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 : 187 : & 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 : 187 : & 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 : 187 : write(std_out,'(a,5x,2a)') ch10,"New overlap for ikpt = ",adjustl(tag)
1541 [ + + ]: 1714 : do itot=1,dimoverlap
1542 : 1714 : write(std_out,'(100f7.3)') (largeoverlap(itot,itot1),itot1=1,dimoverlap)
1543 : : end do
1544 : :
1545 : : ! chipsivect -> chipsi
1546 [ + + ]: 2900 : do ib=1,mbandc
1547 : 2419 : itot = 0
1548 [ + + ]: 9105 : do iatom=1,natom
1549 : 6499 : lpawu = paw_dmft%lpawu(iatom)
1550 [ + + ]: 6499 : if (lpawu == -1) cycle
1551 : 5139 : ndim = nspinor * (2*lpawu+1)
1552 [ + + ]: 67554 : paw_dmft%chipsi(1:ndim,ib,jkpt,isppol,iatom) = mat_tmp(itot+1:itot+ndim,ib)
1553 : 3779 : itot = itot + ndim
1554 : : end do ! iatom
1555 : : end do ! ib
1556 : :
1557 : : ! End big loop over isppol
1558 : : end do !isppol
1559 : :
1560 [ - + ]: 107 : ABI_FREE(chipsivect)
1561 [ - + ]: 107 : ABI_FREE(largeoverlap)
1562 [ - + ]: 107 : ABI_FREE(mat_tmp)
1563 : :
1564 : : end if ! option
1565 : :
1566 [ - + - - : 186 : 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 : 100 : 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 : 100 : integer, allocatable :: displs(:),recvcounts(:)
1593 : 100 : complex(dp), allocatable :: buffer(:),buffer_tot(:)
1594 : : !************************************************************************
1595 : :
1596 : 100 : me_kpt = paw_dmft%distrib%me_kpt
1597 : 100 : mkmem = paw_dmft%distrib%nkpt_mem(me_kpt+1)
1598 : 100 : mbandc = paw_dmft%mbandc
1599 : 100 : natom = paw_dmft%natom
1600 : 100 : nproc = paw_dmft%nproc
1601 : 100 : nspinor = paw_dmft%nspinor
1602 : 100 : nsppol = paw_dmft%nsppol
1603 : 100 : shift = paw_dmft%distrib%shiftk
1604 : :
1605 : 100 : siz_buf = 0
1606 [ + + ]: 439 : do iatom=1,natom
1607 : 339 : lpawu = paw_dmft%lpawu(iatom)
1608 [ + + ]: 339 : if (lpawu == -1) cycle
1609 : 439 : siz_buf = siz_buf + 2*lpawu + 1
1610 : : end do ! iatom
1611 : 100 : siz_buf = siz_buf * nspinor * mbandc * nsppol
1612 : :
1613 [ - + + - : 300 : ABI_MALLOC(recvcounts,(nproc))
+ - + - -
+ - + ]
1614 [ - + + - : 200 : ABI_MALLOC(displs,(nproc))
- + - + ]
1615 : :
1616 [ + + ]: 428 : recvcounts(:) = siz_buf * paw_dmft%distrib%nkpt_mem(:)
1617 : 100 : displs(1) = 0
1618 [ + + ]: 328 : do irank=2,nproc
1619 : 328 : displs(irank) = displs(irank-1) + recvcounts(irank-1)
1620 : : end do ! irank
1621 : :
1622 [ + + + - : 300 : ABI_MALLOC(buffer,(recvcounts(me_kpt+1)))
+ - + + -
+ - + ]
1623 [ - + + - : 300 : ABI_MALLOC(buffer_tot,(recvcounts(nproc)+displs(nproc)))
+ - + - -
+ - + ]
1624 : :
1625 : 100 : ibuf = 0
1626 [ + + ]: 428 : do ikpt=1,mkmem
1627 [ + + ]: 1329 : do iatom=1,natom
1628 : 901 : lpawu = paw_dmft%lpawu(iatom)
1629 [ + + ]: 901 : if (lpawu == -1) cycle
1630 : 374 : ndim = nspinor * (2*lpawu+1)
1631 [ + + ]: 1259 : do isppol=1,nsppol
1632 [ + + ]: 10495 : do ib=1,mbandc
1633 [ + + ]: 92078 : buffer(ibuf+1:ibuf+ndim) = paw_dmft%chipsi(1:ndim,ib,ikpt+shift,isppol,iatom)
1634 : 9594 : ibuf = ibuf + ndim
1635 : : end do ! ib
1636 : : end do ! isppol
1637 : : end do ! iatom
1638 : : end do ! ikpt
1639 : :
1640 : 100 : call xmpi_allgatherv(buffer(:),recvcounts(me_kpt+1),buffer_tot(:),recvcounts(:),displs(:),paw_dmft%distrib%comm_kpt,ierr)
1641 : :
1642 : 100 : ibuf = 0
1643 [ + + ]: 917 : do ikpt=1,paw_dmft%nkpt
1644 [ + + ]: 3105 : do iatom=1,natom
1645 : 2188 : lpawu = paw_dmft%lpawu(iatom)
1646 [ + + ]: 2188 : if (lpawu == -1) cycle
1647 : 977 : ndim = nspinor * (2*lpawu+1)
1648 [ + + ]: 3170 : do isppol=1,nsppol
1649 [ + + ]: 31750 : do ib=1,mbandc
1650 [ + + ]: 315092 : paw_dmft%chipsi(1:ndim,ib,ikpt,isppol,iatom) = buffer_tot(ibuf+1:ibuf+ndim)
1651 : 29562 : ibuf = ibuf + ndim
1652 : : end do ! ib
1653 : : end do ! isppol
1654 : : end do ! iatom
1655 : : end do ! ikpt
1656 : :
1657 : 100 : ABI_FREE(recvcounts)
1658 : 100 : ABI_FREE(displs)
1659 [ - + ]: 100 : ABI_FREE(buffer)
1660 [ - + ]: 100 : ABI_FREE(buffer_tot)
1661 : :
1662 : 100 : 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 : 38 : 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 : 38 : type(oper_type) :: ham_a
1698 : 38 : type(oper_type) :: ham_b
1699 : 38 : type(oper_type) :: ham_squarelocal
1700 : 38 : type(oper_type) :: ham_squareks
1701 : : integer :: iband1,iband2,ikpt,isppol
1702 : : !************************************************************************
1703 : :
1704 : : ! call init_oper(paw_dmft,self_minus_hdc)
1705 : 38 : call init_oper(paw_dmft,ham_a)
1706 : 38 : call init_oper(paw_dmft,ham_b)
1707 : 38 : call init_oper(paw_dmft,ham_squareks)
1708 : 38 : 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 [ + + ]: 677 : do iband1=1,paw_dmft%mbandc
1729 [ + + ]: 32996 : do iband2=1,paw_dmft%mbandc
1730 [ + + ]: 534228 : do ikpt=1,paw_dmft%nkpt
1731 [ + + ]: 1044849 : do isppol=1,paw_dmft%nsppol
1732 [ + + ]: 1012530 : if(iband1==iband2) then
1733 : 10980 : ham_a%ks(iband1,iband2,ikpt,isppol) = cmplx(paw_dmft%eigen_dft(iband1,ikpt,isppol),0.d0,kind=dp)
1734 : : else
1735 : 500280 : 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 : 38 : call downfold_oper(ham_a,paw_dmft)
1745 : :
1746 : : ! Symetrise the local quantity (energy levels)
1747 : : !---------------------------------------------
1748 : 38 : 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 : 38 : call copy_oper(ham_a,ham_b)
1753 : :
1754 : : ! Compute ham_squareks%ks : Make a product of the two KS version
1755 : : !------------------------------------------------------------------
1756 : 38 : call prod_oper(ham_a,ham_b,ham_squareks,1)
1757 : :
1758 : : ! Compute ham_squareks%matlu
1759 : : !---------------------------
1760 : 38 : call downfold_oper(ham_squareks,paw_dmft)
1761 : :
1762 : : ! Symetrise ham_squareks%matlu
1763 : : !------------------------------
1764 : 38 : 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 : 38 : call prod_oper(ham_a,ham_b,ham_squarelocal,2)
1774 : :
1775 : : ! Compute the product in local orbitals
1776 : : !--------------------------------------
1777 : 38 : 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 : 38 : 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 : 38 : 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 : 38 : call destroy_oper(ham_squarelocal)
1801 : : ! call destroy_oper(self_minus_hdc)
1802 : 38 : call destroy_oper(ham_a)
1803 : 38 : call destroy_oper(ham_b)
1804 : 38 : call destroy_oper(ham_squareks)
1805 : :
1806 : :
1807 [ - + - - : 38 : 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 : : !!***
|
