TABLE OF CONTENTS


ABINIT/mlwfovlp_proj [ Functions ]

[ Top ] [ Functions ]

NAME

 mlwfovlp_proj

FUNCTION

 Routine which computes projection A_{mn}(k)
 for Wannier code (www.wannier.org f90 version).

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (BAmadon,FJollet,TRangel,drh)
 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

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
                                          and each |p_lmn> non-local projector
  dtset <type(dataset_type)>=all input variables for this dataset
  filew90_win = secondary input file for wannier90   (WAS NOT USED IN v6.7.1 - so has been temporarily removed)
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  lproj= flag 0: no projections, 1: random projections,
              2: projections on atomic orbitals
              3: projections on projectors
  mband=maximum number of bands
  mkmem =number of k points treated by this node.
  npwarr(nkpt)=number of planewaves in basis at this k point
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw.
  natom=number of atoms in cell.
  nattyp(ntypat)= # atoms of each type.
  nkpt=number of k points.
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat=number of types of atoms in unit cell.
  num_bands=number of bands actually used to construct the wannier function
  nwan= number of wannier fonctions (read in wannier90.win).
  proj_l(mband)= angular part of the projection function (quantum number l)
  proj_m(mband)= angular part of the projection function (quantum number m)
  proj_radial(mband)= radial part of the projection.
  proj_site(3,mband)= site of the projection.
  proj_x(3,mband)= x axis for the projection.
  proj_z(3,mband)= z axis for the projection.
  proj_zona(mband)= extension of the radial part.
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down

OUTPUT

  A_matrix(num_bands,nwan,nkpt,nsppol)= Matrix of projections needed by wannier_run
  ( also wannier90random.amn is written)

SIDE EFFECTS

  (only writing, printing)

NOTES

PARENTS

      mlwfovlp

CHILDREN

      mlwfovlp_radial,mlwfovlp_ylmfac,wrtout,ylm_cmplx

SOURCE

 66 #if defined HAVE_CONFIG_H
 67 #include "config.h"
 68 #endif
 69 
 70 #include "abi_common.h"
 71 
 72  subroutine mlwfovlp_proj(A_matrix,band_in,cg,cprj,dtset,gprimd,just_augmentation,kg,&
 73 &lproj,max_num_bands,mband,mkmem,mpi_enreg,mpw,mwan,natom,nattyp,&
 74 &nkpt,npwarr,nspinor,&
 75 &nsppol,ntypat,num_bands,nwan,pawtab,proj_l,proj_m,proj_radial,&
 76 &proj_site,proj_x,proj_z,proj_zona,psps,spin,ucvol)
 77 
 78  use defs_basis
 79  use defs_datatypes
 80  use defs_abitypes
 81  use defs_wannier90
 82  use m_errors
 83  use m_profiling_abi
 84  use m_xmpi
 85  use m_sort
 86 
 87  use m_numeric_tools, only : uniformrandom
 88  use m_pawtab,  only : pawtab_type
 89  use m_pawcprj, only : pawcprj_type
 90  use m_paw_sphharm, only : ylm_cmplx
 91 
 92 !This section has been created automatically by the script Abilint (TD).
 93 !Do not modify the following lines by hand.
 94 #undef ABI_FUNC
 95 #define ABI_FUNC 'mlwfovlp_proj'
 96  use interfaces_14_hidewrite
 97  use interfaces_67_common, except_this_one => mlwfovlp_proj
 98 !End of the abilint section
 99 
100  implicit none
101 
102 !Arguments ------------------------------------
103 !scalars
104  complex(dpc),parameter :: c1=(1._dp,0._dp)
105  integer,intent(in) :: lproj,max_num_bands,mband,mkmem,mpw,mwan,natom,nkpt,nspinor,nsppol
106  integer,intent(in) :: ntypat,spin
107  type(MPI_type),intent(in) :: mpi_enreg
108  type(dataset_type),intent(in) :: dtset
109  type(pseudopotential_type),intent(in) :: psps
110 !arrays
111  integer ::nattyp(ntypat)
112  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt),num_bands(nsppol),nwan(nsppol),proj_l(mband,nsppol)
113  integer,intent(in) :: proj_m(mband,nsppol)
114  integer,intent(inout)::proj_radial(mband,nsppol)
115  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
116  real(dp),intent(in) :: gprimd(3,3),proj_site(3,mband,nsppol)
117  real(dp),intent(in) :: proj_x(3,mband,nsppol),proj_z(3,mband,nsppol),proj_zona(mband,nsppol)
118  complex(dpc),intent(out) :: A_matrix(max_num_bands,mwan,nkpt,nsppol)
119 !character(len=fnlen),intent(in) :: filew90_win(nsppol)
120  logical,intent(in) :: band_in(mband,nsppol)
121  logical,intent(in)::just_augmentation(mwan,nsppol)
122  type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
123  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
124 
125 !Local variables-------------------------------
126 !scalars
127  integer :: iatom,iatprjn,iband,iband1,iband2,ibg,icat,icg,icg_shift
128  integer :: idum,idx,ikg,ikpt,ilmn,ipw,iproj
129  integer :: ispinor,isppol,itypat,iwan,jband,jj1,libprjn
130  integer :: lmn_size,natprjn,nband_k,nbprjn,npw_k
131  integer :: sumtmp
132  integer :: max_lmax,max_lmax2,mproj,nprocs,spaceComm,rank
133  real(dp),parameter :: qtol=2.0d-8
134  real(dp) :: arg,norm_error,norm_error_bar
135  real(dp) :: ucvol,x1,x2,xnorm,xnormb,xx,yy,zz
136  complex(dpc) :: amn_tmp(nspinor)
137  complex(dpc) :: cstr_fact
138  character(len=500) :: message
139 !arrays
140  integer :: kg_k(3,mpw),lmax(nsppol),lmax2(nsppol),nproj(nsppol)
141  integer,allocatable :: lprjn(:),npprjn(:)
142  real(dp) :: kpg(3),kpt(3)
143  real(dp),allocatable :: amn(:,:,:,:,:),amn2(:,:,:,:,:,:,:)
144  real(dp),allocatable :: gsum2(:),kpg2(:),radial(:)
145  complex(dpc),allocatable :: gf(:,:),gft_lm(:)
146  complex(dpc),allocatable :: ylmc_fac(:,:,:),ylmcp(:)
147 
148 !no_abirules
149 !Tables 3.1 & 3.2, User guide
150  integer,save :: orb_l_defs(-5:3)=(/2,2,1,1,1,0,1,2,3/)
151 ! integer,parameter :: mtransfo(0:3,7)=&
152 !&  reshape((/1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,-2,-1,2,1,0,0,0,-1,1,2,-2,-3,3/),(/4,7/))
153 
154 !************************************************************************
155 
156 
157 !mpi initialization
158  spaceComm=MPI_enreg%comm_cell
159  nprocs=xmpi_comm_size(spaceComm)
160  rank=MPI_enreg%me_kpt
161 
162 !Check input variables
163  if ((lproj/=1).and.(lproj/=2).and.(lproj/=5)) then
164    write(message, '(3a)' )&
165 &   ' Value of lproj no allowed ',ch10,&
166 &   ' Action : change lproj.'
167    MSG_ERROR(message)
168  end if
169 
170  write(message, '(a,a)' )ch10,&
171 & '** mlwfovlp_proj:  compute A_matrix of initial guess for wannier functions'
172  call wrtout(std_out,message,'COLL')
173 
174 !Initialize to 0.d0
175  A_matrix(:,:,:,:)=cmplx(0.d0,0.d0)
176 
177 
178 !
179 !End of preliminarities
180 !
181 
182 !********************* Write Random projectors
183  if(lproj==1) then
184    idum=123456
185 !  Compute random projections
186    ABI_ALLOCATE(amn,(2,mband,mwan,nkpt,nsppol))
187    amn=zero
188    do isppol=1,nsppol
189      if(spin.ne.0 .and. spin.ne.isppol) cycle
190      do ikpt=1,nkpt
191 !
192 !      MPI: cycle over kpts not treated by this node
193 !
194        if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
195 !      write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') ikpt,rank
196 
197 !
198        do iband1=1,mband
199          xnormb=0.d0
200          do iband2=1,nwan(isppol)
201            x1=uniformrandom(idum)
202            x2=uniformrandom(idum)
203            xnorm=sqrt(x1**2+x2**2)
204            xnormb=xnormb+xnorm
205            amn(1,iband1,iband2,ikpt,isppol)=x1
206            amn(2,iband1,iband2,ikpt,isppol)=x2
207          end do
208          do iband2=1,nwan(isppol)
209            amn(1,iband1,iband2,ikpt,isppol)=amn(1,iband1,iband2,ikpt,isppol)/xnormb
210            amn(2,iband1,iband2,ikpt,isppol)=amn(2,iband1,iband2,ikpt,isppol)/xnormb
211          end do !iband2
212        end do !iband1
213      end do !ikpt
214    end do !isppol
215    do isppol=1,nsppol
216      if(spin.ne.0 .and. spin.ne.isppol) cycle
217      do ikpt=1,nkpt
218 !
219 !      MPI: cycle over kpts not treated by this node
220 !
221        if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
222 !
223        do iband2=1,nwan(isppol)
224          jband=0
225          do iband1=1,mband
226            if(band_in(iband1,isppol)) then
227              jband=jband+1
228              if(jband.gt.num_bands(isppol)) then
229                write(message, '(3a)' )&
230 &               '  Value of jband is above num_bands ',ch10,&
231 &               '  Action : contact Abinit group'
232                MSG_ERROR(message)
233              end if
234              A_matrix(jband,iband2,ikpt,isppol)=cmplx(amn(1,iband1,iband2,ikpt,isppol),amn(2,iband1,iband2,ikpt,isppol))
235            end if
236          end do !iband1
237        end do !iband2
238      end do !ikpt
239    end do !isppol
240    ABI_DEALLOCATE(amn)
241  end if
242 
243 !********************* Projection on atomic orbitals based on .win file
244  if( lproj==2) then !based on .win file
245    nproj(:)=nwan(:)/nspinor !if spinors, then the number of projections are
246    mproj=maxval(nproj(:))
247 !  half the total of wannier functions
248 !
249 !  obtain lmax and lmax2
250    lmax(:)=0
251    lmax2(:)=0
252 !
253    do isppol=1,nsppol
254      if(spin.ne.0 .and. spin.ne.isppol) cycle
255      do iproj=1,nproj(isppol)
256        lmax(isppol)=max(lmax(isppol),orb_l_defs(proj_l(iproj,isppol)))
257      end do !iproj
258      lmax2(isppol)=(lmax(isppol)+1)**2
259    end do !isppol
260    max_lmax=maxval(lmax(:))
261    max_lmax2=maxval(lmax2(:))
262 !  Allocate arrays
263    ABI_ALLOCATE(ylmc_fac,(max_lmax2,mproj,nsppol))
264 !
265 !  get ylmfac, factor used for rotations and hybrid orbitals
266    do isppol=1,nsppol
267      if(spin.ne.0 .and. spin.ne.isppol) cycle
268      call mlwfovlp_ylmfac(ylmc_fac(1:lmax2(isppol),1:nproj(isppol),isppol),lmax(isppol),lmax2(isppol),&
269 &     mband,nproj(isppol),proj_l(:,isppol),proj_m(:,isppol),proj_x(:,:,isppol),proj_z(:,:,isppol))
270    end do
271 !
272    norm_error=zero
273    norm_error_bar=zero
274    icg=0
275 !
276    do isppol=1,nsppol
277 !    Allocate arrays
278      if(spin.eq.0 .or. spin.eq.isppol) then
279 !      this has to be done this way because the variable icg changes at the end of the
280 !      cycle. We cannot just skip the hole cycle.
281        ABI_ALLOCATE(gf,(mpw,nproj(isppol)))
282        ABI_ALLOCATE(gft_lm,(lmax2(isppol)))
283        ABI_ALLOCATE(gsum2,(nproj(isppol)))
284        ABI_ALLOCATE(kpg2,(mpw))
285        ABI_ALLOCATE(radial,(lmax2(isppol)))
286        ABI_ALLOCATE(ylmcp,(lmax2(isppol)))
287      end if
288 !
289      ikg=0
290      do ikpt=1, nkpt
291 !
292 !      MPI: cycle over kpts not treated by this node
293 !
294        if (ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-rank)/=0) CYCLE
295 !
296        if(spin.eq.0 .or. spin.eq.isppol) then
297          write(message, '(a,i6,a,2i6)' ) &
298 &         '   processor',rank,' will compute k-point,spin=',ikpt,isppol
299          call wrtout(std_out,  message,'COLL')
300        end if
301 !
302 !      Initialize variables
303        npw_k=npwarr(ikpt)
304        gsum2(:)=0.d0
305        gf(:,:) = (0.d0,0.d0)
306        kpt(:)=dtset%kpt(:,ikpt)
307        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
308 
309        do ipw=1, npw_k
310          kpg(1)= (kpt(1) + real(kg_k(1,ipw),dp))     !k+G
311          kpg(2)= (kpt(2) + real(kg_k(2,ipw),dp))
312          kpg(3)= (kpt(3) + real(kg_k(3,ipw),dp))
313 !
314 !        Calculate modulus of k+G
315          xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3)
316          yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3)
317          zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3)
318          kpg2(ipw)= two_pi*sqrt(xx**2+yy**2+zz**2)
319 !
320 !        Complex Y_lm for k+G
321          if(lmax(isppol)==0) then
322            ylmcp(1)=c1/sqrt(four_pi)
323          else
324            call ylm_cmplx(lmax(isppol),ylmcp,xx,yy,zz)
325          end if
326 !
327          if(spin.eq.0 .or. spin.eq.isppol) then
328 !          !
329            do iproj=1,nproj(isppol)
330 !
331 !            In PAW, we can use proj_radial > 4 to indicate that we just
332 !            want the in-sphere contribution
333 !
334              if( psps%usepaw==1) then
335                if( just_augmentation(iproj,isppol)) cycle
336              end if
337 !
338 !            obtain radial part
339              call mlwfovlp_radial(proj_zona(iproj,isppol),lmax(isppol),lmax2(isppol)&
340 &             ,radial,proj_radial(iproj,isppol),kpg2(ipw))
341 !
342 !            scale complex representation of projector orbital with radial functions
343 !            of appropriate l
344              gft_lm(:)=radial(:)*ylmc_fac(1:lmax2(isppol),iproj,isppol)
345 !
346 !            complex structure factor for projector orbital position
347              arg = ( kpg(1)*proj_site(1,iproj,isppol) + &
348 &             kpg(2)*proj_site(2,iproj,isppol) + &
349 &             kpg(3)*proj_site(3,iproj,isppol) ) * 2*pi
350              cstr_fact = cmplx(cos(arg), -sin(arg) )
351 !
352 !            obtain guiding functions
353              gf(ipw,iproj)=cstr_fact*dot_product(ylmcp,gft_lm)
354 !
355              gsum2(iproj)=gsum2(iproj)+real(gf(ipw,iproj))**2+aimag(gf(ipw,iproj))**2
356            end do !iproj
357          end if !spin
358        end do !ipw
359 !
360        if(spin.eq.0 .or. spin.eq.isppol) then
361          do iproj=1,nproj(isppol)
362 !
363 !          In PAW, we can use proj_radial > 4 to indicate that we just
364 !          want the in-sphere contribution
365 !
366            if(psps%usepaw==1 ) then
367              if (just_augmentation(iproj,isppol)) cycle
368            end if
369 !
370            gsum2(iproj)=16._dp*pi**2*gsum2(iproj)/ucvol
371            gf(:,iproj)=gf(:,iproj)/sqrt(gsum2(iproj))
372            norm_error=max(abs(gsum2(iproj)-one),norm_error)
373            norm_error_bar=norm_error_bar+(gsum2(iproj)-one)**2
374          end do !iproj
375 !
376 !        Guiding functions are computed.
377 !        compute overlaps of gaussian projectors and wave functions
378          do iproj=1,nproj(isppol)
379 !
380 !          In PAW, we can use proj_radial > 4 to indicate that we just
381 !          want the in-sphere contribution
382 !
383            if(psps%usepaw==1 ) then
384              if ( just_augmentation(iproj,isppol)) cycle
385            end if
386 !
387            jband=0
388            do iband=1,mband
389              if(band_in(iband,isppol)) then
390                icg_shift=npw_k*nspinor*(iband-1)+icg
391                jband=jband+1
392                amn_tmp(:)=cmplx(0.d0,0.d0)
393                do ispinor=1,nspinor
394                  do ipw=1,npw_k
395 !
396 !                  The case of spinors is tricky, we have nproj =  nwan/2
397 !                  so we project to spin up and spin down separately, to have at
398 !                  the end an amn matrix with nwan projections.
399 !
400 !                  idx=ipw*nspinor - (nspinor-ispinor)
401                    idx=ipw+(ispinor-1)*npw_k
402                    amn_tmp(ispinor)=amn_tmp(ispinor)+gf(ipw,iproj)*cmplx(cg(1,idx+icg_shift),-cg(2,idx+icg_shift))
403                  end do !ipw
404                end do !ispinor
405                do ispinor=1,nspinor
406                  iwan=(iproj*nspinor)- (nspinor-ispinor)
407                  A_matrix(jband,iwan,ikpt,isppol)=amn_tmp(ispinor)
408                end do
409              end if !band_in
410            end do !iband
411          end do !iproj
412        end if !spin==isppol
413        icg=icg+npw_k*nspinor*mband
414        ikg=ikg+npw_k
415      end do !ikpt
416 !    Deallocations
417      if(spin.eq.0 .or. spin.eq.isppol) then
418        ABI_DEALLOCATE(gf)
419        ABI_DEALLOCATE(gft_lm)
420        ABI_DEALLOCATE(gsum2)
421        ABI_DEALLOCATE(kpg2)
422        ABI_DEALLOCATE(radial)
423        ABI_DEALLOCATE(ylmcp)
424      end if
425    end do !isppol
426 !
427 !  if(isppol==1) then
428 !  norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)),dp))
429 !  else
430 !  if(spin==0)    norm_error_bar=sqrt(norm_error_bar/real(nkpt*(nwan(1)+nwan(2)),dp))
431 !  if(spin==1)    norm_error_bar=sqrt(norm_error_bar/real(nkpt*nwan(1),dp))
432 !  if(spin==2)    norm_error_bar=sqrt(norm_error_bar/real(nkpt*nwan(2),dp))
433 !  end if
434 !  if(norm_error>0.05_dp) then
435 !  write(message, '(6a,f6.3,a,f6.3,12a)' )ch10,&
436 !  &     ' mlwfovlp_proj : WARNING',ch10,&
437 !  &     '  normalization error for wannier projectors',ch10,&
438 !  &     '  is',norm_error_bar,' (average) and',norm_error,' (max).',ch10,&
439 !  &     '  this may indicate more cell-to-cell overlap of the radial functions',ch10,&
440 !  &     '  than you want.',ch10,&
441 !  &     '  Action : modify zona (inverse range of radial functions)',ch10,&
442 !  '  under "begin projectors" in ',trim(filew90_win),' file',ch10
443 !  call wrtout(std_out,message,'COLL')
444 !  end if
445 !
446 !  !Deallocate
447 !  deallocate(ylmc_fac)
448 !
449    ABI_DEALLOCATE(ylmc_fac)
450  end if !lproj==2
451 
452 
453 !*************** computes projection  from PROJECTORS ********************
454  if(lproj==3) then  !! if LPROJPRJ
455 !  ----- set values for projections --------------------- ! INPUT
456 !  nbprjn:number of  different l-values for projectors
457 !  lprjn: value of l for each projectors par ordre croissant
458 !  npprjn: number of projectors for each lprjn
459    natprjn=1  ! atoms with wannier functions are first
460    if(natprjn/=1) then ! in this case lprjn should depend on iatprjn
461      MSG_ERROR("natprjn/=1")
462    end if
463    nbprjn=2
464    ABI_ALLOCATE(lprjn,(nbprjn))
465    lprjn(1)=0
466    lprjn(2)=1
467    ABI_ALLOCATE(npprjn,(0:lprjn(nbprjn)))
468    npprjn(0)=1
469    npprjn(1)=1
470 !  --- test coherence of nbprjn and nwan
471    sumtmp=0
472    do iatprjn=1,natprjn
473      do libprjn=0,lprjn(nbprjn)
474        sumtmp=sumtmp+(2*libprjn+1)*npprjn(libprjn)
475      end do
476    end do
477    if(sumtmp/=nwan(1)) then
478      write(std_out,*) "Number of Wannier orbitals is not equal to number of projections"
479      write(std_out,*) "Action: check values of lprjn,npprjn % nwan"
480      write(std_out,*) "nwan, sumtmp=",nwan,sumtmp
481      MSG_ERROR("Aborting now")
482    end if
483 !  --- end test of coherence
484    ABI_ALLOCATE(amn2,(2,natom,nsppol,nkpt,mband,nspinor,nwan(1)))
485    if(psps%usepaw==1) then
486      amn2=zero
487      ibg=0
488      do isppol=1,nsppol
489        do ikpt=1,nkpt
490          nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
491          do iband=1,nband_k
492 !          write(std_out,*)"amn2",iband,ibg,ikpt
493            do ispinor=1,nspinor
494              icat=1
495              do itypat=1,dtset%ntypat
496                lmn_size=pawtab(itypat)%lmn_size
497                do iatom=icat,icat+nattyp(itypat)-1
498                  jj1=0
499                  do ilmn=1,lmn_size
500                    if(iatom.le.natprjn) then
501 !                    do iwan=1,nwan
502                      do libprjn=0,lprjn(nbprjn)
503 !                      if (psps%indlmn(1,ilmn,itypat)==proj_l(iwan)) then
504 !                      if (psps%indlmn(2,ilmn,itypat)==mtransfo(proj_l(iwan),proj_m(iwan))) then
505                        if (psps%indlmn(1,ilmn,itypat)==libprjn) then
506                          if (psps%indlmn(3,ilmn,itypat)<=npprjn(libprjn)) then
507                            if(band_in(iband,isppol)) then
508                              jj1=jj1+1
509                              if(jj1>nwan(isppol)) then
510                                write(std_out,*) "number of wannier orbitals is lower than lmn_size"
511                                write(std_out,*) jj1,nwan(isppol)
512                                MSG_ERROR("Aborting now")
513                              end if
514                              amn2(1,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(1,ilmn)
515                              amn2(2,iatom,isppol,ikpt,iband,ispinor,jj1)=cprj(iatom,iband+ibg)%cp(2,ilmn)
516                            end if
517                          end if
518                        end if
519                      end do ! libprjn
520 !                    endif
521 !                    endif
522 !                    enddo ! iwan
523                    end if ! natprjn
524                  end do !ilmn
525                end do ! iatom
526                icat=icat+nattyp(itypat)
527              end do ! itypat
528            end do ! ispinor
529          end do !iband
530          ibg=ibg+nband_k*nspinor
531 !        write(std_out,*)'amn2b',iband,ibg,ikpt
532        end do !ikpt
533      end do ! isppol
534 
535 !    -----------------------  Save Amn   --------------------
536      do isppol=1,nsppol
537        do ikpt=1,nkpt
538          do iband2=1,nwan(isppol)
539            jband=0
540            do iband1=1,mband
541              if(band_in(iband1,isppol)) then
542                jband=jband+1
543                A_matrix(jband,iband2,ikpt,isppol)=&
544 &               cmplx(amn2(1,1,1,ikpt,iband1,1,iband2),amn2(2,1,1,ikpt,iband1,1,iband2))
545              end if
546            end do
547          end do
548        end do
549      end do
550    end if !usepaw
551    ABI_DEALLOCATE(amn2)
552    ABI_DEALLOCATE(npprjn)
553    ABI_DEALLOCATE(lprjn)
554 
555  end if ! lproj==3
556 
557 
558 end subroutine mlwfovlp_proj