TABLE OF CONTENTS


ABINIT/datafordmft [ Functions ]

[ Top ] [ Functions ]

NAME

 datafordmft

FUNCTION

  Compute psichi (and print some data for check)

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors .

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
        -gprimd(3,3)=dimensional reciprocal space primitive translations
        -indsym(4,nsym,natom)=indirect indexing array for atom labels
        -symrec(3,3,nsym)=symmetry operations in reciprocal space
        - nsym= number of symetry operations
  cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
                                          and each |p_lmn> non-local projector
  dimcprj(natom) = dimension for cprj
  dtset <type(dataset_type)>=all input variables for this dataset
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  fermie= Fermi energy
  lda_occup <type(oper_type)> = occupations in the correlated orbitals in LDA
  mband=maximum number of bands
  mkmem =number of k points treated by this node
  mpi_enreg=information about MPI parallelization
  nkpt=number of k points.
  my_nspinor=number of spinorial components of the wavefunctions (on current proc)
  nsppol=1 for unpolarized, 2 for spin-polarized
  occ(mband*nkpt*nsppol) = occupancies of KS states.
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang)>=paw angular mesh and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  unpaw = file number for cprj

OUTPUT

  paw_dmft%psichi(nsppol,nkpt,mband,my_nspinor,dtset%natom,(2*maxlpawu+1))): projections <Psi|chi>
  paw_dmft%eigen(paw_dmft%nsppol,paw_dmft%nkpt,paw_dmft%mband)

SIDE EFFECTS

  (only writing, printing)

NOTES

PARENTS

      outscfcv,vtorho

CHILDREN

SOURCE

104 subroutine datafordmft(cryst_struc,cprj,dimcprj,dtset,eigen,fermie,&
105 & lda_occup,mband,mband_cprj,mkmem,mpi_enreg,nkpt,my_nspinor,nsppol,occ,&
106 & paw_dmft,paw_ij,pawang,pawtab,psps,usecprj,unpaw,nbandkss)
107 
108  use defs_basis
109  use defs_datatypes
110  use defs_abitypes
111  use defs_wvltypes
112  use m_abicore
113  use m_errors
114  use m_xmpi
115 
116  use m_io_tools,  only : open_file
117  use m_matlu, only: matlu_type,init_matlu,sym_matlu,copy_matlu,print_matlu,diff_matlu,destroy_matlu
118  use m_crystal, only : crystal_t
119  use m_oper, only : oper_type
120 
121  use m_pawang, only : pawang_type
122  use m_pawtab, only : pawtab_type
123  use m_paw_ij, only : paw_ij_type
124  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free
125  use m_paw_dmft, only: paw_dmft_type
126  use m_mpinfo,   only : proc_distrb_cycle
127 
128 !This section has been created automatically by the script Abilint (TD).
129 !Do not modify the following lines by hand.
130 #undef ABI_FUNC
131 #define ABI_FUNC 'datafordmft'
132 !End of the abilint section
133 
134  implicit none
135 
136 !Arguments ------------------------------------
137 !scalars
138  integer,intent(in) :: mband,mband_cprj,mkmem
139  integer,intent(in) :: nkpt,my_nspinor,nsppol
140  integer,intent(in) :: unpaw,usecprj
141  integer, optional, intent(in) :: nbandkss
142  real(dp),intent(in) :: fermie
143  type(MPI_type),intent(in) :: mpi_enreg
144  type(dataset_type),intent(in) :: dtset
145  type(oper_type), intent(inout) :: lda_occup !vz_i
146  type(pawang_type),intent(in) :: pawang
147  type(pseudopotential_type),intent(in) :: psps
148  type(crystal_t),intent(in) :: cryst_struc
149 !arrays
150  integer, intent(in) :: dimcprj(cryst_struc%natom)
151  real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
152  real(dp),intent(in) :: occ(mband*nkpt*nsppol)
153  type(paw_ij_type),intent(in) :: paw_ij(:)
154 ! type(pawcprj_type) :: cprj(cryst_struc%natom,my_nspinor*mband*mkmem*nsppol)
155  type(pawcprj_type), intent(in) :: cprj(cryst_struc%natom,my_nspinor*mband_cprj*mkmem*nsppol*usecprj)
156  type(paw_dmft_type), intent(inout) :: paw_dmft
157  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
158 !Local variables-------------------------------
159 !scalars
160  integer :: band_index,dimpsichi,facpara
161  integer :: iat,iatom,ib,iband,ibandc,ibg,icat,icount_proj_ilmn,idijeff,ierr,ierrr,ikpt
162  integer :: ilmn,im,im1,iorder_cprj,ispinor,ispinor1,isppol,itypat,ilmn1
163  integer :: jj1,ldim,lmn_size,lpawu
164  integer :: m1,m1_t2g,m1_t2g_mod,maxnproju,me,natom,nband_k,nband_k_cprj
165  integer :: nbandi,nbandf,nnn,nprocband,nsploop,option,opt_renorm,spaceComm,unt
166  real(dp) :: ph0phiint_used
167  character(len=500) :: message
168 !arrays
169  real(dp) :: chinorm
170  complex(dpc), allocatable :: buffer1(:)
171  type(matlu_type), allocatable :: loc_occ_check(:)
172  type(matlu_type), allocatable :: loc_norm_check(:)
173  type(matlu_type), allocatable :: xocc_check(:)
174  type(matlu_type), allocatable :: xnorm_check(:)
175  type(matlu_type), allocatable :: matlu_temp(:)
176  logical :: lprojchi,t2g
177  integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
178  type(pawcprj_type),allocatable :: cwaveprj(:,:)
179 
180 !************************************************************************
181 
182 !DBG_ENTER("COLL")
183 !Fake test to keep fermie as argument. REMOVE IT AS SOON AS POSSIBLE ...
184  if(fermie>huge(zero))chinorm=zero
185 
186  facpara=1 !mpi_enreg%nproc
187  if(abs(dtset%pawprtvol)>=3) then
188    write(message,*) " number of k-points used is nkpt=nkpt ", nkpt
189    call wrtout(std_out,  message,'COLL')
190    write(message,*) " warning: parallelised version        ", nkpt
191    call wrtout(std_out,  message,'COLL')
192    write(message,*) " weights k-points used is wtk=wtk"
193    call wrtout(std_out,  message,'COLL')
194  end if
195 
196  if(usecprj==0) then
197    write(message,*) "  usecprj=0 : BUG in datafordmft",usecprj
198    MSG_BUG(message)
199  end if
200 
201  if(my_nspinor/=dtset%nspinor) then
202    write(message,*) "  my_nspinor=/dtset%nspinor, datafordmft not working in this case",&
203 &   my_nspinor,dtset%nspinor
204    MSG_ERROR(message)
205  end if
206 
207 
208 !do ib=1,my_nspinor*mband_cprj*mkmem*nsppol*usecprj
209 !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
210 !enddo
211 
212 !----------------------------------- MPI-------------------------------------
213 
214 !Init parallelism
215  spaceComm=mpi_enreg%comm_cell
216  if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
217  me=mpi_enreg%me_kpt
218  iorder_cprj=0
219  ABI_CHECK(dtset%mkmem/=0,"mkmem==0 not supported anymore!")
220 !todo_ab: extract cprj from file unpaw in the following..
221 !call abi_abort('COLL')
222 
223 !----------------------------------- MPI-------------------------------------
224 
225  nbandi=paw_dmft%dmftbandi
226  nbandf=paw_dmft%dmftbandf
227  lprojchi=.false.
228  lprojchi=.true.
229  t2g=(paw_dmft%dmftqmc_t2g==1)
230  natom=cryst_struc%natom
231 
232 !if(mpi_enreg%me==0) write(7886,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc
233 !if(mpi_enreg%me==1) write(7887,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc
234 !if(mpi_enreg%me==2) write(7888,*) "in datafordmft", mpi_enreg%me, mpi_enreg%nproc
235  write(message,'(2a)') ch10,&
236 & '  == Prepare data for DMFT calculation  '
237  call wrtout(std_out,message,'COLL')
238  if(abs(dtset%pawprtvol)>=3) then
239    write(message, '(a,a)' ) ch10,&
240 &   '---------------------------------------------------------------'
241 !  call wrtout(ab_out,message,'COLL')
242 !  call wrtout(std_out,  message,'COLL')
243    call wrtout(std_out,  message,'COLL')
244    write(message, '(a,a,a,a,a,a,a,a,a,a,a,a)' ) ch10,&
245 &   '  Print useful data (as a check)',ch10,&
246 &   '  - Overlap of KS wfc with atomic orbital inside sphere',ch10,&
247 &   '  - Eigenvalues',ch10,&
248 &   '  - Weights of k-points',ch10,&
249 &   '  - Number of spins ',ch10,&
250 &   '  - Number of states'
251 !  call wrtout(ab_out,message,'COLL')
252 !  call wrtout(std_out,  message,'COLL')
253    call wrtout(std_out,  message,'COLL')
254    write(message, '(a,a)' ) ch10,&
255 &   '---------------------------------------------------------------'
256  end if
257  if(dtset%nstep==0.and.dtset%nbandkss==0) then
258    message = 'nstep should be greater than 1'
259    MSG_BUG(message)
260  end if
261 
262 !********************* Max Values for U terms.
263 !maxlpawu=0
264  maxnproju=0
265  do iatom=1,natom
266    if(pawtab(dtset%typat(iatom))%lpawu.ne.-1 .and. pawtab(dtset%typat(iatom))%nproju.gt.maxnproju)&
267 &   maxnproju=pawtab(dtset%typat(iatom))%nproju
268  end do
269 !*****************   in forlb.eig
270  if(me.eq.0.and.abs(dtset%pawprtvol)>=3) then
271    if (open_file('forlb.eig',message,newunit=unt,form='formatted',status='unknown') /= 0) then
272      MSG_ERROR(message)
273    end if
274    rewind(unt)
275    write(unt,*) "Number of bands,   spins, and  k-point; and spin-orbit flag"
276    write(unt,*) mband,nsppol,nkpt,my_nspinor,nbandi,nbandf
277    write(unt,*) " For each k-point, eigenvalues for each band"
278    write(unt,*) (dtset%wtk(ikpt)*facpara,ikpt=1,nkpt)
279    band_index=0
280    do isppol=1,nsppol
281      write(unt,*) " For spin"
282      write(unt,*)  isppol
283      do ikpt=1,nkpt
284        nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
285        ibandc=0
286        write(unt,*) " For k-point"
287        write(unt,*)  ikpt
288        do iband=1,mband
289          if(paw_dmft%band_in(iband)) then
290            ibandc=ibandc+1
291            write(unt, '(2i6,4x,f20.15)' ) ibandc,ikpt,eigen(iband+band_index)*2.d0
292          end if
293        end do
294        band_index=band_index+nband_k
295      end do
296    end do
297    close(unt)
298  end if ! proc=me
299 
300 !==   put eigen into eigen_lda
301  band_index=0
302  do isppol=1,nsppol
303    do ikpt=1,nkpt
304      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
305      ibandc=0
306      do iband=1,mband
307        if(paw_dmft%band_in(iband)) then
308          ibandc=ibandc+1
309          paw_dmft%eigen_lda(isppol,ikpt,ibandc)=eigen(iband+band_index) ! in Ha
310         ! paw_dmft%eigen_lda(isppol,ikpt,ibandc)=fermie
311        end if
312      end do
313      band_index=band_index+nband_k
314    end do
315  end do
316 
317  if(abs(dtset%pawprtvol)>=3) then
318    write(message, '(a,a)' ) ch10,&
319 &   '   datafordmft :  eigenvalues written'
320    call wrtout(std_out,  message,'COLL')
321  end if
322 !==========================================================================
323 !***************** Compute  <Psi|Chi>=\sum_{proja} <Psi|P_a><phi_a|Chi>
324 !==========================================================================
325 !write(std_out,*) "size(cprj,dim=1)",size(cprj,dim=1),size(cprj,dim=2),dtset%mband,dtset%mkmem,dtset%nkpt
326 
327 !Allocate temporary cwaveprj storage
328  ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,my_nspinor))
329 !write(std_out,*) "before alloc cprj"
330 !write(std_out,*) size(cwaveprj,dim=1),size(cwaveprj,dim=2),size(dimcprj,dim=1)
331 
332  call pawcprj_alloc(cwaveprj,0,dimcprj)
333 !write(std_out,*) "after alloc cprj"
334 
335  nprocband=(mband/mband_cprj)
336  ibg=0
337  do isppol=1,nsppol
338    do ikpt=1,nkpt
339      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
340      nband_k_cprj=nband_k/nprocband
341 
342      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
343 
344 !    write(2011,*) ikpt
345 !    ibsp=ibg
346      ibandc=0
347 !    LOOP OVER BANDS
348      ib=0
349      do iband=1,nband_k
350 
351        if(paw_dmft%band_in(iband)) ibandc=ibandc+1
352 
353 !      Parallelization: treat only some bands
354        if (dtset%paral_kgb==1) then
355          if (mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band)/=mpi_enreg%me_band) cycle
356        else
357          if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle
358        end if
359        ib=ib+1
360 
361        do ispinor=1,my_nspinor
362 !        ibsp =~ (ikpt-1)*nband*my_nspinor+iband
363 !        ibsp=ibsp+1
364          icat=1
365 !        write(std_out,*) isppol,ikpt,iband,ispinor
366          iat=0 ! to declare
367          do itypat=1,dtset%ntypat
368            lmn_size=pawtab(itypat)%lmn_size
369 !          write(std_out,*) isppol,ikpt,iband,ispinor
370            do iatom=icat,icat+cryst_struc%nattyp(itypat)-1
371              lpawu=pawtab(dtset%typat(iatom))%lpawu
372 !            ----------   t2g case
373              if(paw_dmft%dmftqmc_t2g==1.and.lpawu/=-1) then
374                if(lpawu==2) then
375 !                lpawu==2 must be chosen in input and thus in
376 !                pawtab. On the contrary, paw_dmft now has
377 !                lpawu=1
378                  m1_t2g=0 ! index for psichi which has a dimension 3
379                else
380                  write(message,'(a,a,i4,i4,2a)')  ch10,&
381 &                 '  Wrong use of dmftqmc_t2g',paw_dmft%dmftqmc_t2g,lpawu,ch10,&
382 &                 ' Action: desactivate qmftqmc_t2g or use lpawu=1'
383                  MSG_ERROR(message)
384                end if
385              end if
386 !            ----------   t2g case
387 !            if(isppol==2) write(std_out,*) "ee",size(cprj(iatom,ibsp)%cp(:,:))
388 
389              iat=iat+1
390              jj1=0
391              if(lpawu.ne.-1) then
392 
393                call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,ib,ibg,ikpt,&
394 &               iorder_cprj,isppol,mband_cprj,mkmem,natom,1,nband_k_cprj,&
395 &               my_nspinor,nsppol,unpaw,mpicomm=mpi_enreg%comm_kpt,&
396 &               proc_distrb=mpi_enreg%proc_distrb)
397 !              write(std_out,*) "M-2",isppol,ikpt,iband,iatom,&
398 !              &             (cwaveprj(iatom,ispinor)%cp(1,13)**2+cwaveprj(iatom,ispinor)%cp(2,13)**2) ,ibsp
399 
400 !              chinorm=(pawtab(itypat)%phiphjint(1))
401                chinorm=1.d0
402 !              write(std_out,*) isppol,ikpt,iband,ispinor,iat
403                do ilmn=1,lmn_size
404 !                write(std_out,*) ilmn
405 !                ------------ Select l=lpawu.
406                  if (psps%indlmn(1,ilmn,itypat)==lpawu) then
407 !                  ------------ Check that the band is choosen (within nbandi and nbandf)
408                    if(paw_dmft%band_in(iband)) then
409 !                    if(ilmn==13) then
410 !                    write(std_out,*) "M-2c",isppol,ikpt,iband,iatom
411 !                    write(std_out,*) "M-2b",isppol,ikpt,ibandc,&
412 !                    &             (cwaveprj(iatom,ispinor)%cp(1,13)**2+cwaveprj(iatom,ispinor)%cp(2,13)**2)
413 !                    endif
414 !                    write(std_out,*) "inside paw_dmft%band_in",iband
415                      jj1=jj1+1
416                      if(jj1>pawtab(dtset%typat(iatom))%nproju*(2*lpawu+1)) then
417                        write(message,'(a,a,a,a)')  ch10,&
418 &                       ' jj1 is not correct in datafordmft',ch10,&
419 &                       ' Action: CONTACT Abinit group'
420                        MSG_BUG(message)
421                      end if ! jj1
422                      icount_proj_ilmn=psps%indlmn(3,ilmn,itypat)  ! n: nb of projector
423 !                    write(std_out,*) "icount_proj_ilmn",icount_proj_ilmn
424                      m1=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1
425 !                    ---- if lprochi=true do the sum over every projectors
426 !                    ---- if lprochi=false: . do the sum over only ground state projector
427 !                    ----                   . and take ph0phiint(icount_proj_ilmn)=1
428 
429 !                    call pawcprj_get(cryst_struc%atindx1,cwaveprj,cprj,natom,&
430 !                    &             ib,ibg,ikpt,iorder_cprj,isppol,mband_cprj,mkmem,&
431 !                    &             natom,1,nband_k_cprj,my_nspinor,nsppol,unpaw,&
432 !                    &             mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
433 
434                      if(lprojchi.or.icount_proj_ilmn==1) then
435                        if(.not.lprojchi) ph0phiint_used=one
436                        if(lprojchi) ph0phiint_used=pawtab(itypat)%ph0phiint(icount_proj_ilmn)
437                        if(paw_dmft%dmftqmc_t2g==1) then ! t2g case
438 
439                          if(m1==1.or.m1==2.or.m1==4) then ! t2g case
440                            m1_t2g=m1_t2g+1  ! t2g case1
441                            m1_t2g_mod=mod(m1_t2g-1,3)+1
442 !                          write(std_out,*) "M0",isppol,ikpt,iband,ilmn,cprj(iatom,ibsp)%cp(1,ilmn)
443 !                          write(std_out,*) "M1",m1,m1_t2g,iband,ilmn,icount_proj_ilmn
444                            paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g_mod)=&  ! t2g case
445 &                          paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g_mod)+&  ! t2g case
446 !                          &                           cmplx(cprj(iatom,ibsp)%cp(1,ilmn)*ph0phiint_used,&  ! t2g case
447 !                          &                           cprj(iatom,ibsp)%cp(2,ilmn)*ph0phiint_used,kind=dp)  ! t2g case
448 &                          cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,&  ! t2g case
449 &                          cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp)  ! t2g case
450                          end if  !t2g case
451                        else
452                          paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)=&
453 &                         paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)+&
454 !                        &                         cmplx(cprj(iatom,ibsp)%cp(1,ilmn)*ph0phiint_used,&
455 !                        &                         cprj(iatom,ibsp)%cp(2,ilmn)*ph0phiint_used,kind=dp)
456 &                         cmplx(cwaveprj(iatom,ispinor)%cp(1,ilmn)*ph0phiint_used,&
457 &                         cwaveprj(iatom,ispinor)%cp(2,ilmn)*ph0phiint_used,kind=dp)
458 
459 !                        if(ibandc==3.and.iat==1.and.m1==1) then
460 !                        write(std_out,'(a,3i5)') "psichi integers",isppol,ikpt,ispinor
461 !                        write(std_out,'(a,2i5,2e16.7)') "psichi IB3 iAT1 IM1",ilmn,icount_proj_ilmn,&
462 !                        &             real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1))
463 !                        write(std_out,'(a,2i5,2e16.7)') "cwaveprj IB3 iAT1 IM1",ilmn,icount_proj_ilmn,cwaveprj(iatom,ispinor)%cp(1,ilmn) &
464 !                        &                    , cwaveprj(iatom,ispinor)%cp(2,ilmn)
465 !                        endif
466                        end if
467                      end if ! lprojchi=.true. (always)
468                    end if ! nband belongs to paw_dmft%band_in
469                  end if ! L=lpawu
470                end do !ilmn : loop over L,M,N
471              end if ! If lpawu/=-1
472            end do ! iatom
473            icat=icat+cryst_struc%nattyp(itypat)
474          end do ! itypat
475        end do ! ispinor
476      end do !iband
477      ibg=ibg+nband_k_cprj*my_nspinor
478 !    bdtot_index=bdtot_index+nband_k ! useless  (only for occ)
479    end do !ikpt
480  end do ! isppol
481 !do isppol=1,nsppol
482 !do ikpt=1,nkpt
483 !do ispinor=1,my_nspinor
484 !write(std_out,*) "psichi integers",isppol,ikpt,ispinor
485 !write(std_out,*) "psichi IB3 iAT1 IM1",&
486 !&             real(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1)), imag(paw_dmft%psichi(isppol,ikpt,3,ispinor,1,1))
487 !
488 !enddo
489 !enddo
490 !enddo
491 !call abi_abort('COLL')
492  if(abs(dtset%pawprtvol)>=3) then
493    write(message,*) "chinorm used here =",chinorm
494    call wrtout(std_out,  message,'COLL')
495  end if
496 
497 !deallocate temporary cwaveprj/cprj storage
498  call pawcprj_free(cwaveprj)
499  ABI_DATATYPE_DEALLOCATE(cwaveprj)
500 
501 !==========================================================================
502 !********************* Gather information for MPI before printing
503 !==========================================================================
504 
505  dimpsichi=2*nsppol*nkpt*mband*my_nspinor*natom*(2*paw_dmft%maxlpawu+1)
506  ABI_ALLOCATE(buffer1,(dimpsichi))
507  nnn=0
508 !write(176,*) "beg",psichi
509  do isppol=1,nsppol
510    do ikpt=1,nkpt
511      do ibandc=1,paw_dmft%mbandc
512        do ispinor=1,my_nspinor
513          do iat=1,natom
514            do m1=1,2*paw_dmft%maxlpawu+1
515 !            do m=1,2
516              nnn=nnn+1
517              buffer1(nnn)=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)
518 !            enddo
519            end do
520          end do
521        end do
522      end do
523    end do
524  end do
525  call xmpi_barrier(spaceComm)
526  call xmpi_sum(buffer1,spaceComm ,ierr)
527  if (dtset%paral_kgb==1.and.nprocband>1) then
528    call xmpi_sum(buffer1,mpi_enreg%comm_band,ierr) !Build sum over band processors
529  end if
530  call xmpi_barrier(spaceComm)
531  nnn=0
532  do isppol=1,nsppol
533    do ikpt=1,nkpt
534      do ibandc=1,paw_dmft%mbandc
535        do ispinor=1,my_nspinor
536          do iat=1,natom
537            do m1=1,2*paw_dmft%maxlpawu+1
538 !            do m=1,2
539              nnn=nnn+1
540              paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1)=buffer1(nnn)
541 !            enddo
542            end do
543          end do
544        end do
545      end do
546    end do
547  end do
548  ABI_DEALLOCATE(buffer1)
549 !write(177,*) "end",psichi
550 
551 !do isppol=1,nsppol
552 !do ikpt=1,nkpt
553 !do ibandc=1,paw_dmft%mbandc
554 !do ispinor=1,my_nspinor
555 !write(std_out,*) "psichigather",isppol,ikpt,ibandc,&
556 !&             real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2+&
557 !&             imag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,1,1))**2
558 !
559 !enddo
560 !enddo
561 !enddo
562 !enddo
563 
564  call xmpi_barrier(spaceComm)
565 !if(mpi_enreg%me.eq.0) write(177,*) "end",psichi
566 !if(mpi_enreg%me.eq.1) write(178,*) "end",psichi
567 !if(mpi_enreg%me.eq.2) write(179,*) "end",psichi
568  call xmpi_barrier(spaceComm)
569 !==========================================================================
570 !********* WRITE psichi in file for reference
571 !==========================================================================
572  if(me.eq.0) then
573    call psichi_print(dtset,cryst_struc%nattyp,cryst_struc%ntypat,nkpt,my_nspinor,&
574 &   nsppol,paw_dmft,pawtab,psps,t2g)
575  end if ! proc=me
576 !==========================================================================
577 !********************* Check normalization  and occupations ***************
578 !==========================================================================
579  ABI_DATATYPE_ALLOCATE(xocc_check,(natom))
580  ABI_DATATYPE_ALLOCATE(xnorm_check,(natom))
581  call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,xocc_check)
582  call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,xnorm_check)
583  call psichi_check(dtset,cryst_struc%nattyp,nkpt,my_nspinor,&
584 & nsppol,cryst_struc%ntypat,paw_dmft,pawtab,psps,xocc_check,xnorm_check)
585 !==========================================================================
586 !***************  write checks  *******************************************
587 !==========================================================================
588  if(abs(dtset%pawprtvol)>=3) then
589    write(message,*) "normalization computed"
590    call wrtout(std_out,  message,'COLL')
591  end if
592 
593  ABI_DATATYPE_ALLOCATE(loc_occ_check,(natom))
594  ABI_DATATYPE_ALLOCATE(loc_norm_check,(natom))
595  call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,loc_occ_check)
596  call init_matlu(natom,my_nspinor,nsppol,paw_dmft%lpawu,loc_norm_check)
597  call copy_matlu(xocc_check,loc_occ_check,natom)
598  call copy_matlu(xnorm_check,loc_norm_check,natom)
599 
600  write(message,'(2a,i4)')  ch10," == Check: Occupations and Norm from psichi are"
601  call wrtout(std_out,  message,'COLL')
602 
603  if(paw_dmft%dmftcheck>=1) then
604 !  print occupations
605    write(message,'(2a,i4)')  ch10,'  ------ Unsymetrised Occupation'
606    call wrtout(std_out,  message,'COLL')
607 
608    call print_matlu(xocc_check,natom,dtset%pawprtvol)
609 
610 !  print norms
611    write(message,'(2a,i4)')  ch10,'  ------ Unsymetrised Norm'
612    call wrtout(std_out,  message,'COLL')
613 
614    call print_matlu(xnorm_check,natom,dtset%pawprtvol)
615  end if
616 
617 !symetrise and print occupations
618  call sym_matlu(cryst_struc,loc_occ_check,pawang)
619 
620  write(message,'(2a,i4)')  ch10,'  ------ Symetrised Occupation'
621  call wrtout(std_out,  message,'COLL')
622 
623  call print_matlu(loc_occ_check,natom,dtset%pawprtvol)
624 
625 !symetrise and print norms
626  call sym_matlu(cryst_struc,loc_norm_check,pawang)
627 
628  write(message,'(2a,i4)')  ch10,'  ------ Symetrised Norm'
629  call wrtout(std_out,  message,'COLL')
630 
631  call print_matlu(loc_norm_check,natom,dtset%pawprtvol)
632 
633 !deallocations
634  do iatom=1,natom
635    lda_occup%matlu(iatom)%mat=loc_occ_check(iatom)%mat
636  end do
637 
638 !Tests density matrix LDA+U and density matrix computed here.
639  if(paw_dmft%dmftcheck==2.or.(paw_dmft%dmftbandi==1)) then
640    ABI_DATATYPE_ALLOCATE(matlu_temp,(natom))
641    call init_matlu(natom,paw_dmft%nspinor,paw_dmft%nsppol,paw_dmft%lpawu,matlu_temp)
642    do iatom=1,natom
643      if(paw_dmft%lpawu(iatom).ne.-1) then
644        ldim=2*paw_dmft%lpawu(iatom)+1
645        nsploop=max(paw_dmft%nsppol,paw_dmft%nspinor**2)
646        do idijeff=1,nsploop
647          ispinor=0
648          ispinor1=0
649          if(nsploop==2) then
650            isppol=spinor_idxs(1,idijeff)
651            ispinor=1
652            ispinor1=1
653          else if(nsploop==4) then
654            isppol=1
655            ispinor=spinor_idxs(1,idijeff)
656            ispinor1=spinor_idxs(2,idijeff)
657          else if(nsploop==1) then
658            isppol=1
659            ispinor=1
660            ispinor1=1
661          else
662            write(message,'(2a)') " BUG in datafordmft: nsploop should be equal to 2 or 4"
663            call wrtout(std_out,message,'COLL')
664          end if
665          do im1 = 1 , ldim
666            do im = 1 ,  ldim
667              if(my_nspinor==2) matlu_temp(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=&
668 &             cmplx(paw_ij(iatom)%noccmmp(1,im,im1,idijeff),paw_ij(iatom)%noccmmp(2,im,im1,idijeff),kind=dp)
669              if(my_nspinor==1) matlu_temp(iatom)%mat(im,im1,isppol,ispinor,ispinor1)=&
670 &             cmplx(paw_ij(iatom)%noccmmp(1,im,im1,idijeff),zero,kind=dp)
671            end do
672          end do
673        end do
674      end if
675    end do
676    if(paw_dmft%dmftcheck==2) option=1
677    if(paw_dmft%dmftcheck<=1) option=0
678    call diff_matlu("LDA+U density matrix from INPUT wfk",&
679 &   "Direct calculation of density matrix with psichi from DIAGONALIZED wfk",&
680 &   matlu_temp,lda_occup%matlu,natom,option,tol3,ierrr) !tol1 tol2 tol3
681    if(ierrr==-1) then
682      write(message,'(10a)') ch10,&
683 &     '    -> These two quantities should agree if three conditions are fulfilled',ch10,&
684 &     '         -  input wavefunctions come from the same Hamiltonien (e.g LDA/GGA)',ch10,&
685 &     '         -  dmatpuopt is equal to 1',ch10,&
686 &     '         -  all valence states are in the valence',ch10,&
687 &     '    (for experts users: it is not compulsary that these conditions are fulfilled)'
688      call wrtout(std_out,message,'COLL')
689    end if
690 !  write(message,'(2a)') ch10,&
691 !  &   '  ***** => Calculations of density matrices with projections and in LDA+U are coherent****'
692 !  call wrtout(std_out,message,'COLL')
693 
694    call destroy_matlu(matlu_temp,natom)
695    ABI_DATATYPE_DEALLOCATE(matlu_temp)
696  else
697    write(message,'(2a)') ch10,&
698 &   '  Warning: Consistency of density matrices computed from projection has not been checked: use dmftcheck>=2 '
699    call wrtout(std_out,message,'COLL')
700  end if
701 
702  call destroy_matlu(loc_norm_check,natom)
703  ABI_DATATYPE_DEALLOCATE(loc_norm_check)
704  call destroy_matlu(loc_occ_check,natom)
705  ABI_DATATYPE_DEALLOCATE(loc_occ_check)
706 
707  call destroy_matlu(xocc_check,natom)
708  call destroy_matlu(xnorm_check,natom)
709  ABI_DATATYPE_DEALLOCATE(xocc_check)
710  ABI_DATATYPE_DEALLOCATE(xnorm_check)
711 
712  if(present(nbandkss)) then
713    if(me.eq.0.and.nbandkss/=0) then
714 !     opt_renorm=1 ! if ucrpa==1, no need for individual orthonormalization
715      opt_renorm=3
716      if(dtset%ucrpa>=1) opt_renorm=2
717      call psichi_renormalization(cryst_struc,paw_dmft,pawang,opt=opt_renorm)
718      call psichi_print(dtset,cryst_struc%nattyp,cryst_struc%ntypat,nkpt,my_nspinor,&
719 &     nsppol,paw_dmft,pawtab,psps,t2g)
720    end if ! proc=me
721  end if

ABINIT/m_datafordmft [ Modules ]

[ Top ] [ Modules ]

NAME

  m_datafordmft

FUNCTION

 This module produces inputs for the DMFT calculation

COPYRIGHT

 Copyright (C) 2006-2018 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

25 #if defined HAVE_CONFIG_H
26 #include "config.h"
27 #endif
28 
29 #include "abi_common.h"
30 
31 MODULE m_datafordmft
32 
33  use defs_basis
34 
35  implicit none
36 
37  private
38 
39  public :: datafordmft
40  public :: compute_levels
41  public :: psichi_renormalization
42  public :: hybridization_asymptotic_coefficient

m_datafordmft/compute_levels [ Functions ]

[ Top ] [ m_datafordmft ] [ Functions ]

NAME

 compute_levels

FUNCTION

 Compute levels for ctqmc

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

OUTPUT

NOTES

PARENTS

      m_hubbard_one,qmc_prep_ctqmc

CHILDREN

      checkdiag_matlu,loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

1060  subroutine compute_levels(cryst_struc,energy_level,hdc,pawang,paw_dmft,nondiag)
1061 
1062  use defs_basis
1063  use defs_datatypes
1064  use defs_abitypes
1065  use m_errors
1066  use m_abicore
1067 
1068  use m_pawang, only : pawang_type
1069  use m_crystal, only : crystal_t
1070  use m_paw_dmft, only : paw_dmft_type
1071  use m_oper, only : oper_type,loc_oper
1072  use m_matlu, only : sym_matlu, print_matlu, checkdiag_matlu
1073 
1074 !This section has been created automatically by the script Abilint (TD).
1075 !Do not modify the following lines by hand.
1076 #undef ABI_FUNC
1077 #define ABI_FUNC 'compute_levels'
1078 !End of the abilint section
1079 
1080  implicit none
1081 
1082 !Arguments ------------------------------------
1083 !scalars
1084  type(crystal_t),intent(in) :: cryst_struc
1085  type(pawang_type), intent(in) :: pawang
1086  type(oper_type), intent(in) :: hdc
1087  type(paw_dmft_type), intent(in)  :: paw_dmft
1088  type(oper_type),intent(inout)  :: energy_level !vz_i
1089  logical, optional, intent(out) :: nondiag
1090 
1091 !Local variables ------------------------------
1092 ! scalars
1093  integer :: iatom,iband,ikpt,im1,isppol,ispinor
1094  integer :: lpawu,mbandc,natom,nspinor,nsppol,nkpt
1095  character(len=500) :: message
1096 ! arrays
1097 !************************************************************************
1098 
1099  mbandc=paw_dmft%mbandc
1100  nkpt=paw_dmft%nkpt
1101  nsppol=paw_dmft%nsppol
1102  nspinor=paw_dmft%nspinor
1103  natom=paw_dmft%natom
1104  if(present(nondiag)) nondiag=.false.
1105 
1106 !========================
1107 !Get KS eigenvalues
1108 !========================
1109  do iband=1,mbandc
1110    do ikpt=1,nkpt
1111      do isppol=1,nsppol
1112 !      Take \epsilon_{nks}
1113 !      ========================
1114        energy_level%ks(isppol,ikpt,iband,iband)=paw_dmft%eigen_lda(isppol,ikpt,iband)
1115      end do
1116    end do
1117  end do
1118 
1119 
1120 !======================================================================
1121 !Compute atomic levels from projection of \epsilon_{nks} and symetrize
1122 !======================================================================
1123  call loc_oper(energy_level,paw_dmft,1)
1124 ! write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels before sym and only LDA"
1125 ! call wrtout(std_out,message,'COLL')
1126 ! call print_matlu(energy_level%matlu,natom,1)
1127  do iatom = 1 , natom
1128    lpawu=paw_dmft%lpawu(iatom)
1129    if(lpawu/=-1) then
1130      do isppol=1,nsppol
1131        do ispinor=1,nspinor
1132          do im1=1,2*lpawu+1
1133            energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)= &
1134 &           energy_level%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)&
1135 &           -hdc%matlu(iatom)%mat(im1,im1,isppol,ispinor,ispinor)-paw_dmft%fermie
1136          end do
1137        end do
1138      end do
1139 !    write(std_out,*) "DC,fermie",hdc%matlu(iatom)%mat(1,1,1,1,1),paw_dmft%fermie
1140    end if
1141  end do ! natom
1142  call sym_matlu(cryst_struc,energy_level%matlu,pawang)
1143  if(present(nondiag)) call checkdiag_matlu(energy_level%matlu,natom,tol7,nondiag)
1144 
1145  write(message,'(a,2x,a,f13.5)') ch10," == Print Energy levels for Fermi Level=",paw_dmft%fermie
1146  call wrtout(std_out,message,'COLL')
1147 !call print_oper(energy_level,1,paw_dmft,1)
1148  call print_matlu(energy_level%matlu,natom,1)
1149 
1150 
1151 
1152  end subroutine compute_levels

m_datafordmft/hybridization_asymptotic_coefficient [ Functions ]

[ Top ] [ m_datafordmft ] [ Functions ]

NAME

 hybridization_asymptotic_coefficient

FUNCTION

 Compute some components for the limit of hybridization

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cryst_struc <type(crystal_t)>=crystal structure data
  lda_occup
  paw_dmft =  data for self-consistent LDA+DMFT calculations.
  pawang <type(pawang)>=paw angular mesh and related data
  pawtab <type(pawtab)>

OUTPUT

  paw_dmft =  data for self-consistent LDA+DMFT calculations.

NOTES

PARENTS

      qmc_prep_ctqmc

CHILDREN

      add_matlu,copy_oper,destroy_oper,init_oper,loc_oper,prod_oper,sym_matlu

SOURCE

1633 subroutine hybridization_asymptotic_coefficient(cryst_struc,paw_dmft,pawang,hybri_coeff)
1634 
1635 
1636  use defs_basis
1637  use defs_abitypes
1638  use m_errors
1639  use m_abicore
1640  use m_crystal, only : crystal_t
1641  use m_oper, only : oper_type,init_oper,upfold_oper,copy_oper,prod_oper,destroy_oper,loc_oper
1642  use m_matlu, only : matlu_type,init_matlu,add_matlu,destroy_matlu,print_matlu,sym_matlu
1643  use m_paw_dmft, only: paw_dmft_type
1644  use m_pawang, only : pawang_type
1645 
1646 !This section has been created automatically by the script Abilint (TD).
1647 !Do not modify the following lines by hand.
1648 #undef ABI_FUNC
1649 #define ABI_FUNC 'hybridization_asymptotic_coefficient'
1650 !End of the abilint section
1651 
1652  implicit none
1653 
1654 !Arguments ------------------------------------
1655 !scalars
1656  type(crystal_t),intent(in) :: cryst_struc
1657  type(paw_dmft_type), intent(in) :: paw_dmft
1658  type(pawang_type), intent(in) :: pawang
1659  type(matlu_type), intent(inout) :: hybri_coeff(paw_dmft%natom)
1660 !Local variables ------------------------------
1661  type(oper_type)  :: ham_a
1662  type(oper_type)  :: ham_b
1663  type(oper_type)  :: ham_squarelocal
1664  type(oper_type)  :: ham_squareks
1665  integer :: iband1,iband2,ikpt,isppol
1666 !************************************************************************
1667 
1668 ! call init_oper(paw_dmft,self_minus_hdc)
1669  call init_oper(paw_dmft,ham_a)
1670  call init_oper(paw_dmft,ham_b)
1671  call init_oper(paw_dmft,ham_squareks)
1672  call init_oper(paw_dmft,ham_squarelocal)
1673 
1674 ! Create self_minus_hdc%matlu = Sigma-Hdc in local orbitals
1675 ! call add_matlu(self%oper(paw_dmft%dmft_nwlo)%matlu,self%hdc%matlu,&
1676 !&             self_minus_hdc%matlu,cryst_struc%natom,-1)
1677 
1678 !   write(message,'(a,2x,a)') ch10,        "  == self_minus_hdc (1)"
1679 !   call wrtout(std_out,message,'COLL')
1680 !   call print_matlu(self_minus_hdc%matlu,paw_dmft%natom,1,opt_exp=1)
1681 
1682 !! Create self_minus_hdc%ks = Upfold Sigma-Hdc
1683 ! call upfold_oper(self_minus_hdc,paw_dmft,1)
1684 ! call loc_oper(self_minus_hdc,paw_dmft,1)
1685 
1686 !   write(message,'(a,2x,a)') ch10,        "  == self_minus_hdc (2)"
1687 !   call wrtout(std_out,message,'COLL')
1688 !   call print_matlu(self_minus_hdc%matlu,paw_dmft%natom,1,opt_exp=1)
1689 
1690 ! Create ham_a%ks = H_ks  in KS basis
1691 !----------------------------------------------------
1692  do iband1=1,paw_dmft%mbandc
1693    do iband2=1,paw_dmft%mbandc
1694      do ikpt=1,paw_dmft%nkpt
1695        do isppol=1,paw_dmft%nsppol
1696          if(iband1==iband2) then
1697            ham_a%ks(isppol,ikpt,iband1,iband2) = cmplx(paw_dmft%eigen_lda(isppol,ikpt,iband1),0.d0,kind=dp)
1698          else
1699            ham_a%ks(isppol,ikpt,iband1,iband2) = czero
1700          end if
1701        end do
1702      end do
1703    end do
1704  end do
1705 
1706 ! Create ham_a%matlu = H_ks in local orbitals
1707 !---------------------------------------------
1708  call loc_oper(ham_a,paw_dmft,1)
1709 
1710 ! Symetrise the local quantity (energy levels)
1711 !---------------------------------------------
1712  call sym_matlu(cryst_struc,ham_a%matlu,pawang)
1713 
1714 ! Create ham_b%ks : Duplicate both ks and local part of ham_a into ham_b
1715 !-----------------------------------------------------------------------
1716  call copy_oper(ham_a,ham_b)
1717 
1718 ! Compute ham_squareks%ks   : Make a product of the two KS version
1719 !------------------------------------------------------------------
1720  call prod_oper(ham_a,ham_b,ham_squareks,1)
1721 
1722 ! Compute ham_squareks%matlu
1723 !---------------------------
1724  call loc_oper(ham_squareks,paw_dmft,1)
1725 
1726 ! Symetrise ham_squareks%matlu
1727 !------------------------------
1728  call sym_matlu(cryst_struc,ham_squareks%matlu,pawang)
1729 
1730 !   write(message,'(a,2x,a)') ch10,        "  == squareks"
1731 !   call wrtout(std_out,message,'COLL')
1732 !   call print_matlu(ham_squareks%matlu,paw_dmft%natom,1,opt_exp=1)
1733 
1734 
1735 ! Compute ham_squarelocal%matlu
1736 !-------------------------------
1737  call prod_oper(ham_a,ham_b,ham_squarelocal,2)
1738 
1739 ! Compute the product in local orbitals
1740 !--------------------------------------
1741  call sym_matlu(cryst_struc,ham_squarelocal%matlu,pawang)
1742 
1743 !   write(message,'(a,2x,a)') ch10,        "  == squarelocal"
1744 !   call wrtout(std_out,message,'COLL')
1745 !   call print_matlu(ham_squarelocal%matlu,paw_dmft%natom,1,opt_exp=1)
1746 
1747 ! The difference of ham_squareks and ham_squarelocal
1748 ! gives the coefficient that we are looking for ( such that F_ij(iw_n) = -C_ij/(iw_n) ).
1749 !----------------------------------------------------------------------------------------
1750  call add_matlu(ham_squareks%matlu,ham_squarelocal%matlu,hybri_coeff,cryst_struc%natom,-1)
1751 
1752  !  write(message,'(a,2x,a)') ch10,        "  == Coeff C_ij before sym"
1753  !  call wrtout(std_out,message,'COLL')
1754  !  call print_matlu(hybri_coeff,paw_dmft%natom,1,opt_exp=1)
1755 
1756 ! Symetrise the local quantity
1757 !------------------------------
1758  call sym_matlu(cryst_struc,hybri_coeff,pawang)
1759 
1760  !  write(message,'(a,2x,a)') ch10,        "  == Coeff C_ij after sym"
1761  !  call wrtout(std_out,message,'COLL')
1762  !  call print_matlu(hybri_coeff,paw_dmft%natom,1,opt_exp=1)
1763 
1764  call destroy_oper(ham_squarelocal)
1765 ! call destroy_oper(self_minus_hdc)
1766  call destroy_oper(ham_a)
1767  call destroy_oper(ham_b)
1768  call destroy_oper(ham_squareks)
1769 
1770 
1771 end subroutine hybridization_asymptotic_coefficient

m_datafordmft/psichi_check [ Functions ]

[ Top ] [ m_datafordmft ] [ Functions ]

NAME

  psichi_check

FUNCTION

  Check psichi: compute occupations

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  maxnproju = maximum number of projector for LDA+U species
  nattyp(ntypat)= # atoms of each type
  mband= number of bands
  nkpt=number of k points
  my_nspinor=number of spinorial components of the wavefunctions (on current proc)
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat= number of species
  paw_dmft <type(paw_dmft)>=paw data for the self-consistency
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials

  OUTPUTS:
  xocc_check(nsppol,my_nspinor,my_nspinor,natom,2*maxlpawu+1,2*maxlpawu+1): density matrix
  xnorm_check(nsppol,my_nspinor,my_nspinor,natom,2*maxlpawu+1,2*maxlpawu+1): matrix of norms

SIDE EFFECTS

  check psichi: compute norm and occupations

PARENTS

      datafordmft

CHILDREN

SOURCE

 924 subroutine psichi_check(dtset,nattyp,nkpt,my_nspinor,&
 925 & nsppol,ntypat,paw_dmft,pawtab,psps,xocc_check,xnorm_check)
 926 
 927  use m_abicore
 928 
 929  use m_matlu, only: matlu_type,init_matlu,sym_matlu
 930 
 931 !This section has been created automatically by the script Abilint (TD).
 932 !Do not modify the following lines by hand.
 933 #undef ABI_FUNC
 934 #define ABI_FUNC 'psichi_check'
 935 !End of the abilint section
 936 
 937  implicit none
 938 
 939 !Arguments ------------------------------------
 940 !scalars
 941  integer,intent(in) :: nkpt,my_nspinor,nsppol,ntypat
 942 !arrays
 943  integer, intent(in) :: nattyp(ntypat)
 944  type(dataset_type),intent(in) :: dtset
 945  type(pseudopotential_type),intent(in) :: psps
 946  type(paw_dmft_type), intent(in) :: paw_dmft
 947  type(matlu_type), intent(inout) :: xocc_check(dtset%natom) !vz_i
 948  type(matlu_type), intent(inout) :: xnorm_check(dtset%natom) !vz_i
 949  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
 950 !Local variables ------------------------------------
 951  integer :: band_index,facnsppol,ibg,isppol,ikpt,iband,ibandc,ispinor,icat,itypat
 952  integer :: iat,iatom,ilmn,lmn_size,m,m1,nband_k
 953  complex(dpc) :: psichic,psichic1
 954  real(dp) :: chinorm
 955 
 956 ! *********************************************************************
 957    facnsppol=1
 958    if(my_nspinor==1.and.nsppol==1) then
 959      facnsppol=2
 960    end if
 961 
 962    ibg=0
 963    band_index=0
 964    do isppol=1,nsppol
 965      do ikpt=1,nkpt
 966        nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
 967        ibandc=0
 968        do iband=1,nband_k
 969          if(paw_dmft%band_in(iband)) ibandc=ibandc+1
 970          do ispinor=1,my_nspinor
 971            icat=1
 972            iat=0
 973            do itypat=1,dtset%ntypat
 974              lmn_size=pawtab(itypat)%lmn_size
 975              do iatom=icat,icat+nattyp(itypat)-1
 976                iat=iat+1
 977 !              ------------ Select correlated atoms
 978                if(paw_dmft%lpawu(iatom).ne.-1) then
 979                  chinorm=1.d0
 980                  do ilmn=1,lmn_size
 981 !                  ------------ Select l=lpawu.
 982                    if (psps%indlmn(1,ilmn,itypat)==paw_dmft%lpawu(iatom).and.&
 983 &                   psps%indlmn(3,ilmn,itypat)==1) then
 984                      do ilmn1=1,lmn_size
 985 !                      ------------ Select l=lpawu and do not sum over projectors
 986 !                      (this is already done in paw_dmft%psichi)
 987                        if (psps%indlmn(1,ilmn1,itypat)==paw_dmft%lpawu(iatom).and.&
 988 &                       psps%indlmn(3,ilmn1,itypat)==1) then
 989 !                        ------------ Check that the band is choosen (within nbandi and nbandf)
 990                          if(paw_dmft%band_in(iband)) then
 991                            m=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1
 992                            m1=psps%indlmn(2,ilmn1,itypat)+psps%indlmn(1,ilmn,itypat)+1
 993                            if(psps%indlmn(3,ilmn,itypat)==1) then
 994                              do ispinor1=1,my_nspinor
 995                                psichic=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m)
 996                                psichic1=paw_dmft%psichi(isppol,ikpt,ibandc,ispinor1,iat,m1)
 997 !                              ------------ Compute occupation matrix
 998                                xocc_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)=&
 999 &                               xocc_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)&
1000 !                              &               +conjg(psichic)*psichic1*dtset%wtk(ikpt)*facpara*occ(iband+band_index)
1001 &                               +conjg(psichic1)*psichic*dtset%wtk(ikpt)*facpara*occ(iband+band_index)/facnsppol
1002 !                              ------------ Compute norm (should be equal to noccmmp
1003 !                              (dmatpuopt=1) if all bands are taken into account)
1004                                xnorm_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)=&
1005 &                               xnorm_check(iatom)%mat(m,m1,isppol,ispinor,ispinor1)&
1006 !                              &               +conjg(psichic)*psichic1*dtset%wtk(ikpt)*facpara
1007 &                               +conjg(psichic1)*psichic*dtset%wtk(ikpt)*facpara
1008                              end do ! ispinor1
1009                            end if
1010                          end if ! paw_dmft%band_in
1011                        end if
1012                      end do !ilmn1
1013                    end if
1014                  end do !ilmn
1015                end if ! lpawu.ne.-1
1016              end do ! iatom
1017              icat=icat+nattyp(itypat)
1018            end do ! itypat
1019          end do ! ispinor
1020        end do !iband
1021        band_index=band_index+nband_k
1022        ibg=ibg+nband_k*my_nspinor
1023      end do !ikpt
1024    end do ! isppol
1025 
1026  end subroutine psichi_check
1027 !DBG_EXIT("COLL")
1028 end subroutine datafordmft

m_datafordmft/psichi_print [ Functions ]

[ Top ] [ m_datafordmft ] [ Functions ]

NAME

  psichi_print

FUNCTION

  Print psichi for reference

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  maxnproju = maximum number of projector for LDA+U species
  nattyp(ntypat)= # atoms of each type
  mband= number of bands
  nkpt=number of k points
  my_nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  paw_dmft <type(paw_dmft)>=paw data for the self-consistency
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  psichi(2,nsppol,nkpt,mband,my_nspinor,dtset%natom,(2*paw_dmft%maxlpawu+1))) projections for DMFT
  psps <type(pseudopotential_type)>=variables related to pseudopotentials

SIDE EFFECTS

  print psichi in forlb.ovlp

PARENTS

      datafordmft

CHILDREN

SOURCE

757 subroutine psichi_print(dtset,nattyp,ntypat,nkpt,my_nspinor,&
758 &nsppol,paw_dmft,pawtab,psps,t2g)
759 
760  use m_abicore
761  use m_io_tools,  only : open_file
762 
763 !This section has been created automatically by the script Abilint (TD).
764 !Do not modify the following lines by hand.
765 #undef ABI_FUNC
766 #define ABI_FUNC 'psichi_print'
767 !End of the abilint section
768 
769  implicit none
770 !Arguments ------------------------------------
771 !scalars
772  integer,intent(in) :: nkpt,my_nspinor,nsppol,ntypat
773 !arrays
774  integer, intent(in) :: nattyp(ntypat)
775  type(dataset_type),intent(in) :: dtset
776  type(pseudopotential_type),intent(in) :: psps
777  logical t2g
778  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
779  type(paw_dmft_type), intent(in) :: paw_dmft
780 !Local variables ------------------------------------
781  integer :: ibg,isppol,ikpt,iband,ibandc,ispinor,icat,itypat,lmn_size
782  integer :: iat,iatom,jj1,ilmn,m1,nband_k,unt
783  integer :: m1_t2g,ll
784  real(dp) :: chinorm
785  character(len=500) :: msg
786 
787 ! *********************************************************************
788    ll=1
789 
790    if (open_file('forlb.ovlp',msg,newunit=unt,form='formatted',status='unknown') /= 0) then
791      MSG_ERROR(msg)
792    end if
793    rewind(unt)
794 
795 !  Header for calc_uCRPA.F90
796    if  (COUNT(pawtab(:)%lpawu.NE.-1).EQ.1) then
797      do  itypat=1,ntypat
798        if(t2g) then
799          if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",ll,itypat
800        else
801          if(pawtab(itypat)%lpawu.ne.-1) write(unt,*) "l= ",pawtab(itypat)%lpawu,itypat
802        end if
803      end do
804    else
805      write(unt,*) "More than one correlated species"
806    end if
807 
808    write(unt,*) "Bands ",paw_dmft%dmftbandi,paw_dmft%dmftbandf
809 
810    ibg=0
811    do isppol=1,nsppol
812      do ikpt=1,nkpt
813 !      rewind(1023)
814        write(unt,'(a6,2x,i6)') "ikpt =",ikpt
815        nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
816        ibandc=0
817        do iband=1,nband_k
818          if(paw_dmft%band_in(iband)) then
819            ibandc=ibandc+1
820            write(unt,'(a8,2x,i6)') " iband =",iband
821          end if
822          do ispinor=1,my_nspinor
823            icat=1
824 !          write(std_out,*) isppol,ikpt,iband,ispinor
825            iat=0 ! to declare
826            do itypat=1,dtset%ntypat
827              lmn_size=pawtab(itypat)%lmn_size
828 !            write(std_out,*) isppol,ikpt,iband,ispinor
829              do iatom=icat,icat+nattyp(itypat)-1
830                iat=iat+1
831                jj1=0
832                if(pawtab(dtset%typat(iatom))%lpawu.ne.-1) then
833 !                chinorm=(pawtab(itypat)%phiphjint(1))
834                  chinorm=1.d0
835 !                write(std_out,*) isppol,ikpt,iband,ispinor,iat
836                  m1_t2g=0
837                  do ilmn=1,lmn_size
838 !                  write(std_out,*) ilmn
839 !                  ------------ Select l=lpawu.  ---------------------------------------
840                    if (psps%indlmn(1,ilmn,itypat)==pawtab(dtset%typat(iatom))%lpawu.and.psps%indlmn(3,ilmn,itypat)==1) then
841 !                    ------------ Check that the band is choosen (within nbandi and nbandf)
842                      if(paw_dmft%band_in(iband)) then
843                        jj1=jj1+1
844                        if(jj1>(2*pawtab(dtset%typat(iatom))%lpawu+1)) then
845                          write(message,'(a,a,i4,i5,i4)') ch10," Error 2 in datafordmft",jj1,pawtab(dtset%typat(iatom))%lpawu
846                          call wrtout(std_out,  message,'COLL')
847                          MSG_ERROR("Aborting now")
848                        end if ! jj1
849 !                      if(jj1>pawtab(dtset%typat(iatom))%nproju*(2*lpawu+1)) then
850 !                      write(message,'(a,a,a,a)')  ch10,&
851 !                      &                         'BUG: jj1 is not correct in datafordmft psichi_print',ch10,&
852 !                      &                         'Action: CONTACT Abinit group'
853 !                      call wrtout(std_out,  message,'COLL')
854 !                      call abi_abort('COLL')
855 !                      end if ! jj1
856                        m1=psps%indlmn(2,ilmn,itypat)+psps%indlmn(1,ilmn,itypat)+1
857 !                      ----- Print only when the sum over projectors is done
858 !                      write(std_out,*) ilmn,m1
859                        if(t2g) then
860                          if(m1==1.or.m1==2.or.m1==4) then
861                            m1_t2g=m1_t2g+1
862                            write(unt,'(3i6,3x,2f23.15)') isppol, iat, m1,&
863 &                           real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g))/chinorm,&
864 &                           aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1_t2g))/chinorm
865                          end if
866                        else
867                          write(unt,'(3i6,3x,2f23.15)') isppol, iat, m1,&
868 &                         real(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1))/chinorm,&
869 &                         aimag(paw_dmft%psichi(isppol,ikpt,ibandc,ispinor,iat,m1))/chinorm
870                        end if
871 !                      if natom=1 then jj1 maximal value should be 2*lpawu+1
872                      end if ! paw_dmft%band_in
873                    end if
874                  end do !ilmn
875                end if ! lpawu.ne.-1
876              end do ! iatom
877              icat=icat+nattyp(itypat)
878            end do ! itypat
879          end do ! ispinor
880        end do !iband
881        ibg=ibg+nband_k*my_nspinor
882      end do !ikpt
883    end do ! isppol
884 !   write(unt,*) "Fermi level (in Ryd)="
885 !   write(unt,*) fermie*two
886    close(unt)
887  end subroutine psichi_print

m_datafordmft/psichi_renormalization [ Functions ]

[ Top ] [ m_datafordmft ] [ Functions ]

NAME

 psichi_renormalization

FUNCTION

 Renormalize psichi.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (BAmadon)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cryst_struc <type(crystal_t)>= crystal structure data.
  paw_dmft =  data for LDA+DMFT calculations.
  pawang <type(pawang)>=paw angular mesh and related data

OUTPUT

  paw_dmft%psichi(nsppol,nkpt,mband,nspinor,dtset%natom,(2*maxlpawu+1))): projections <Psi|chi> are orthonormalized.

NOTES

PARENTS

      m_datafordmft,m_dmft

CHILDREN

      add_matlu,destroy_oper,gather_matlu,identity_oper,init_oper
      invsqrt_matrix,loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

1189 subroutine psichi_renormalization(cryst_struc,paw_dmft,pawang,opt)
1190 
1191  use defs_basis
1192  use m_errors
1193  use m_abicore
1194 
1195  use m_pawang, only : pawang_type
1196  use m_paw_dmft, only: paw_dmft_type
1197  use m_crystal, only : crystal_t
1198  use m_oper, only : init_oper,oper_type,identity_oper,loc_oper,destroy_oper,diff_oper
1199  use m_matlu, only : matlu_type,sym_matlu,print_matlu
1200 
1201 !This section has been created automatically by the script Abilint (TD).
1202 !Do not modify the following lines by hand.
1203 #undef ABI_FUNC
1204 #define ABI_FUNC 'psichi_renormalization'
1205 !End of the abilint section
1206 
1207  implicit none
1208 
1209 !Arguments ------------------------------------
1210 !scalars
1211  type(crystal_t),intent(in) :: cryst_struc
1212  type(paw_dmft_type), intent(inout)  :: paw_dmft
1213  type(pawang_type), intent(in) :: pawang
1214  integer, optional, intent(in) :: opt
1215 
1216 !Local variables ------------------------------
1217 !scalars
1218  integer :: iatom,ib,ikpt,im,ispinor,isppol,jkpt
1219  integer :: natom,mbandc,ndim,nkpt,nspinor,nsppol,option
1220  real(dp), pointer ::  temp_wtk(:) => null()
1221  real(dp) ::  pawprtvol
1222  type(oper_type) :: norm
1223  type(oper_type) :: oper_temp
1224  character(len=500) :: message
1225 !arrays
1226 ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:)
1227 !************************************************************************
1228 
1229  DBG_ENTER("COLL")
1230 
1231  option=3
1232  if(present(opt)) then
1233    if(opt==2) option=2
1234    if(opt==3) option=3
1235  end if
1236  pawprtvol=2
1237 
1238  nsppol  = paw_dmft%nsppol
1239  nkpt    = paw_dmft%nkpt
1240  mbandc  = paw_dmft%mbandc
1241  natom   = cryst_struc%natom
1242  nspinor = paw_dmft%nspinor
1243 
1244 
1245 !== Normalize psichi
1246  if(option==1) then
1247 !  ====================================
1248 !  == simply renormalize psichi =======
1249 !  ====================================
1250    write(message,'(2a)') ch10," Psichi are renormalized  "
1251    call wrtout(std_out,  message,'COLL')
1252    do isppol=1,nsppol
1253      do ikpt=1,nkpt
1254        do ib=1,mbandc
1255          do iatom=1,natom
1256            if(paw_dmft%lpawu(iatom).ne.-1) then
1257              ndim=2*paw_dmft%lpawu(iatom)+1
1258              do im=1,ndim
1259                do ispinor=1,nspinor
1260 !                write(std_out,*) "psichi1",paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)
1261                  paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)=     &
1262 &                 paw_dmft%psichi(isppol,ikpt,ib,ispinor,iatom,im)/    &
1263 &                 sqrt(real(norm%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)))
1264                end do ! ispinor
1265              end do ! im
1266            end if
1267          end do ! iatom
1268        end do ! ib
1269      end do ! ikpt
1270    end do ! isppol
1271 !  todo_ab introduce correct orthonormalization in the general case.
1272 
1273  else if(option==2) then ! option==2
1274 !  ====================================
1275 !  == renormalize k-point after k-point
1276 !  ====================================
1277    ABI_ALLOCATE(temp_wtk,(1))
1278 
1279    write(message,'(2a,i5)') ch10,' Nb of K-point',nkpt
1280    call wrtout(std_out,message,'COLL')
1281    do jkpt=1,nkpt  ! jkpt
1282      write(message,'(2a,i5)') ch10,'  == Renormalization for K-point',jkpt
1283      call wrtout(std_out,message,'COLL')
1284      temp_wtk(1)=one
1285 
1286      call normalizepsichi(cryst_struc,1,paw_dmft,pawang,temp_wtk,jkpt)
1287    end do ! jkpt
1288    ABI_DEALLOCATE(temp_wtk)
1289 
1290  else if(option==3)  then  ! option==3
1291 !  ====================================
1292 !  == renormalize the sum over k-points
1293 !  ====================================
1294 !  allocate(temp_wtk(nkpt))
1295    temp_wtk=>paw_dmft%wtk
1296    write(message,'(6a)') ch10,'  ====================================== '&
1297 &   ,ch10,'  == Renormalization for all K-points == '&
1298 &   ,ch10,'  ======================================='
1299    call wrtout(std_out,message,'COLL')
1300    call normalizepsichi(cryst_struc,nkpt,paw_dmft,pawang,temp_wtk)
1301 !  deallocate(temp_wtk)
1302 
1303  end if
1304 
1305 !== Change back repr for norm
1306 
1307 !===============================================
1308 !==  Compute norm with new psichi
1309 !===============================================
1310  call init_oper(paw_dmft,norm,nkpt=paw_dmft%nkpt,wtk=paw_dmft%wtk)
1311 !== Build identity for norm%ks (option=1)
1312  call identity_oper(norm,1)
1313 
1314 !== Deduce norm%matlu from norm%ks with new psichi
1315  call loc_oper(norm,paw_dmft,1)
1316 
1317 !== Print norm%matlu unsymetrized with new psichi
1318  if(pawprtvol>2) then
1319    write(message,'(4a,2a)') &
1320 &   ch10,"  == Check: Overlap with renormalized psichi without symetrization is == "
1321    call wrtout(std_out,message,'COLL')
1322    call print_matlu(norm%matlu,natom,prtopt=1)
1323  end if
1324 
1325 
1326 !== Symetrise norm%matlu with new psichi
1327  call sym_matlu(cryst_struc,norm%matlu,pawang)
1328 
1329 !== Print norm%matlu symetrized with new psichi
1330  if(pawprtvol>2) then
1331    write(message,'(4a,2a)') &
1332 &   ch10,"  == Check: Overlap with renormalized psichi and symetrization is =="
1333    call wrtout(std_out,message,'COLL')
1334    call print_matlu(norm%matlu,natom,prtopt=1,opt_diag=-1)
1335  end if
1336 
1337 !== Check that norm is now the identity
1338  call init_oper(paw_dmft,oper_temp)
1339  call identity_oper(oper_temp,2)
1340  call diff_oper('Overlap after renormalization','Identity',&
1341 & norm,oper_temp,1,tol6)
1342  call destroy_oper(oper_temp)
1343 
1344  call destroy_oper(norm)
1345 
1346  paw_dmft%lpsichiortho=1
1347 
1348  DBG_EXIT("COLL")
1349 
1350  CONTAINS
1351 !===========================================================

psichi_renormalization/normalizepsichi [ Functions ]

[ Top ] [ psichi_renormalization ] [ Functions ]

NAME

  normalizepsichi

FUNCTION

  normalizepsichi psichi from norm

INPUTS

SIDE EFFECTS

  change psichi: normalizepsichi it

PARENTS

      psichi_renormalization

CHILDREN

      add_matlu,destroy_oper,gather_matlu,identity_oper,init_oper
      invsqrt_matrix,loc_oper,print_matlu,sym_matlu,wrtout

SOURCE

1375 subroutine normalizepsichi(cryst_struc,nkpt,paw_dmft,pawang,temp_wtk,jkpt)
1376 
1377  use m_abicore
1378 
1379  use defs_basis
1380  use m_errors
1381 
1382  use m_paw_dmft, only: paw_dmft_type
1383  use m_crystal, only : crystal_t
1384  use m_oper, only : init_oper,oper_type,identity_oper,loc_oper,destroy_oper
1385  use m_matlu, only : gather_matlu,sym_matlu,print_matlu,add_matlu
1386  use m_matrix, only : invsqrt_matrix
1387 
1388 !This section has been created automatically by the script Abilint (TD).
1389 !Do not modify the following lines by hand.
1390 #undef ABI_FUNC
1391 #define ABI_FUNC 'normalizepsichi'
1392 !End of the abilint section
1393 
1394  implicit none
1395 
1396 !Arguments ------------------------------------
1397  integer,intent(in) :: nkpt
1398  integer,optional,intent(in) :: jkpt
1399  real(dp),pointer :: temp_wtk(:)
1400 !scalars
1401  type(crystal_t),intent(in) :: cryst_struc
1402  type(paw_dmft_type), intent(inout)  :: paw_dmft
1403  type(pawang_type), intent(in) :: pawang
1404 
1405 !Local variables ------------------------------
1406 !scalars
1407  integer :: diag,iatom,ib,ikpt1,im,im1,ispinor,ispinor1,isppol,isppol1,jc,jc1
1408  integer :: tndim
1409  integer :: natom,mbandc,ndim,nspinor,nsppol
1410  real(dp) :: pawprtvol
1411  type(oper_type) :: norm1,norm2,norm3
1412  character(len=500) :: message
1413  complex(dpc),allocatable :: wan(:,:,:),sqrtmatinv(:,:)
1414  type(coeff2c_type), allocatable :: overlap(:)
1415 !arrays
1416 ! real(dp),allocatable :: e0pde(:,:,:),omegame0i(:)
1417 !************************************************************************
1418    nsppol  = paw_dmft%nsppol
1419    mbandc  = paw_dmft%mbandc
1420    natom   = cryst_struc%natom
1421    nspinor = paw_dmft%nspinor
1422    pawprtvol=3
1423    diag=0
1424 
1425    if(nkpt/=1.and.present(jkpt)) then
1426      message = 'BUG in psichi_normalization'
1427      MSG_ERROR(message)
1428    end if
1429 
1430 !  *********************************************************************
1431    call init_oper(paw_dmft,norm1,nkpt=nkpt,wtk=temp_wtk)
1432 
1433 !  == Build identity for norm1%ks (option=1)
1434    call identity_oper(norm1,1)
1435 
1436    if(nkpt==1.and.present(jkpt)) then
1437      call loc_oper(norm1,paw_dmft,1,jkpt=jkpt)
1438    end if
1439    if(.not.present(jkpt)) then
1440      call loc_oper(norm1,paw_dmft,1)
1441    end if
1442    if(nkpt>1) then
1443      call sym_matlu(cryst_struc,norm1%matlu,pawang)
1444    end if
1445 
1446    if(pawprtvol>2) then
1447      write(message,'(2a)') ch10,'  - Print norm with current psichi '
1448      call wrtout(std_out,message,'COLL')
1449      call print_matlu(norm1%matlu,natom,prtopt=1,opt_exp=1)
1450    end if
1451 !  ==-------------------------------------
1452 !  == Start loop on atoms
1453    ABI_DATATYPE_ALLOCATE(overlap,(natom))
1454    do iatom=1,natom
1455      if(paw_dmft%lpawu(iatom).ne.-1) then
1456        ndim=2*paw_dmft%lpawu(iatom)+1
1457        tndim=nsppol*nspinor*ndim
1458        ABI_ALLOCATE(overlap(iatom)%value,(tndim,tndim))
1459        overlap(iatom)%value=czero
1460      end if
1461    end do
1462 !  ==-------------------------------------
1463 
1464 !  built large overlap matrix
1465    write(message,'(2a)') ch10,'  - Overlap (before orthonormalization) -'
1466    call wrtout(std_out,message,'COLL')
1467    call gather_matlu(norm1%matlu,overlap,cryst_struc%natom,option=1,prtopt=1)
1468    call destroy_oper(norm1)
1469 
1470 
1471 
1472    do iatom=1,natom
1473      if(paw_dmft%lpawu(iatom).ne.-1) then
1474        ndim=2*paw_dmft%lpawu(iatom)+1
1475        tndim=nsppol*nspinor*ndim
1476        ABI_ALLOCATE(sqrtmatinv,(tndim,tndim))
1477 
1478 !      == Compute Inverse Square root of overlap : O^{-0.5}
1479        !!write(message,'(a,1x,a,e21.14,a,e21.14,a)') "overlap", &
1480        !!"(",real(overlap(1)%value),",",aimag(overlap(1)%value),")"
1481        !!call wrtout(std_out,message,'COLL')
1482        if(diag==0) then
1483          call invsqrt_matrix(overlap(iatom)%value,tndim)
1484          sqrtmatinv=overlap(iatom)%value
1485        else
1486          sqrtmatinv(:,:)=czero
1487          do ib=1,tndim
1488            sqrtmatinv(ib,ib)=cone/(sqrt(overlap(iatom)%value(ib,ib)))
1489          end do
1490        end if
1491 
1492 !      == Apply O^{-0.5} on psichi
1493        ABI_ALLOCATE(wan,(nsppol,nspinor,ndim))
1494 !      write(std_out,*) mbandc,nsppol,nspinor,ndim
1495 !      write(std_out,*)  paw_dmft%psichi(1,1,1,1,1,1)
1496        do ikpt=1,nkpt
1497          do ib=1,mbandc
1498            if(present(jkpt)) then
1499              ikpt1=jkpt
1500            else
1501              ikpt1=ikpt
1502            end if
1503            jc=0
1504            wan=czero
1505            do isppol=1,nsppol
1506              do ispinor=1,nspinor
1507                do im=1,ndim
1508 !                write(std_out,*) "psichi", paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)
1509                  jc=jc+1
1510                  jc1=0
1511                  do isppol1=1,nsppol
1512                    do ispinor1=1,nspinor
1513                      do im1=1,ndim
1514                        jc1=jc1+1
1515                        wan(isppol,ispinor,im)= wan(isppol,ispinor,im) &
1516 &                       + paw_dmft%psichi(isppol1,ikpt1,ib,ispinor1,iatom,im1)*sqrtmatinv(jc,jc1)
1517                      end do ! ispinor1
1518                    end do ! isppol1
1519                  end do ! im1
1520                end do ! im
1521              end do ! ispinor
1522            end do !  isppol
1523            do isppol=1,nsppol
1524              do ispinor=1,nspinor
1525                do im=1,ndim
1526                  paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)=wan(isppol,ispinor,im)
1527 !                write(std_out,*) "psichi2", paw_dmft%psichi(isppol,ikpt1,ib,ispinor,iatom,im)
1528                end do ! ispinor
1529              end do ! isppol
1530            end do ! im
1531          end do ! ib
1532        end do ! ikpt
1533        ABI_DEALLOCATE(wan)
1534        ABI_DEALLOCATE(sqrtmatinv)
1535 !      write(std_out,*)  paw_dmft%psichi(1,1,1,1,1,1)
1536 
1537 !      ==-------------------------------------
1538      end if ! lpawu.ne.-1
1539    end do ! iatom
1540 !  == End loop on atoms
1541 !  ==-------------------------------------
1542    do iatom=1,natom
1543      if(paw_dmft%lpawu(iatom).ne.-1) then
1544        ABI_DEALLOCATE(overlap(iatom)%value)
1545      end if
1546    end do
1547    ABI_DATATYPE_DEALLOCATE(overlap)
1548 
1549 !  ======================================================================
1550 !  == Check norm with new psichi.
1551 !  ======================================================================
1552 
1553    call init_oper(paw_dmft,norm1,nkpt=nkpt,wtk=temp_wtk)
1554 
1555    call identity_oper(norm1,1)
1556 
1557    if(nkpt==1.and.present(jkpt)) then
1558      call loc_oper(norm1,paw_dmft,1,jkpt=jkpt)
1559    end if
1560    if(.not.present(jkpt)) then
1561      call loc_oper(norm1,paw_dmft,1)
1562    end if
1563 
1564    if (nkpt>1) then
1565      call sym_matlu(cryst_struc,norm1%matlu,pawang)
1566    end if
1567 
1568    if(pawprtvol>2) then
1569      write(message,'(2a)') ch10,'  - Print norm with new psichi '
1570      call wrtout(std_out,message,'COLL')
1571      call print_matlu(norm1%matlu,natom,prtopt=1)
1572    end if
1573 
1574 !  ======================================================================
1575 !  == Check that norm-identity is zero
1576 !  ======================================================================
1577    call init_oper(paw_dmft,norm2,nkpt=nkpt,wtk=temp_wtk)
1578    call init_oper(paw_dmft,norm3,nkpt=nkpt,wtk=temp_wtk)
1579    call identity_oper(norm2,2)
1580    call add_matlu(norm1%matlu,norm2%matlu,norm3%matlu,natom,-1)
1581    call destroy_oper(norm2)
1582    if(pawprtvol>2) then
1583      write(message,'(2a)') ch10,'  - Print norm with new psichi minus Identity '
1584      call wrtout(std_out,message,'COLL')
1585      call print_matlu(norm3%matlu,natom,prtopt=1,opt_exp=1)
1586    end if
1587    call destroy_oper(norm3)
1588 
1589    call destroy_oper(norm1)
1590 !  call flush(std_out)           ! debug debug  debug   debug
1591 !  MSG_ERROR("Stop for debugging")
1592 
1593  end subroutine normalizepsichi
1594 
1595 end subroutine psichi_renormalization