TABLE OF CONTENTS


ABINIT/mlwfovlp_pw [ Functions ]

[ Top ] [ Functions ]

NAME

 mlwfovlp_pw

FUNCTION

 Routine which computes PW part of overlap M_{mn}(k,b) 
 for Wannier code (www.wannier.org f90 version).

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (BAmadon,FJollet,T Rangel)
 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.
  g1(3,nkpt,nntot) = G vector shift which is necessary to obtain k1+b
  iwav(mband,nkpt,nsppol): shift for pw components in cg.
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mband=maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem =number of k points treated by this node.
  mpi_enreg=informations about MPI parallelization
  mpw=maximum dimensioned size of npw.
  nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
  ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  ovikp(nkpt,nntot)= gives  nntot value of k2 (in the BZ) for each k1  (k2=k1+b mod(G))
  seed_name= seed_name of files containing cg for all k-points to be used with MPI
  spin = just used for nsppol>1 ; 0 both, 1 just spin up, 2 just spin down

OUTPUT

  cm1(2,mband,mband,nntot,nkpt,nsppol): overlap <u_(nk1)|u_(mk1+b)>.

SIDE EFFECTS

  (only writing, printing)

NOTES

PARENTS

      mlwfovlp

CHILDREN

      wrtout,xmpi_barrier,xmpi_sum

SOURCE

 53 #if defined HAVE_CONFIG_H
 54 #include "config.h"
 55 #endif
 56 
 57 #include "abi_common.h"
 58 
 59 subroutine mlwfovlp_pw(cg,cm1,g1,iwav,kg,mband,mkmem,mpi_enreg,mpw,nfft,ngfft,nkpt,nntot,&
 60 &  npwarr,nspinor,nsppol,ovikp,seed_name,spin)
 61 
 62  use defs_basis
 63  use defs_abitypes
 64  use m_profiling_abi
 65  use m_errors
 66  use m_xmpi
 67 
 68 !This section has been created automatically by the script Abilint (TD).
 69 !Do not modify the following lines by hand.
 70 #undef ABI_FUNC
 71 #define ABI_FUNC 'mlwfovlp_pw'
 72  use interfaces_14_hidewrite
 73 !End of the abilint section
 74 
 75  implicit none
 76 
 77 !Arguments ------------------------------------
 78 !scalars
 79  integer,intent(in) :: mband,mkmem,mpw,nfft,nkpt,nntot
 80  integer,intent(in) :: nspinor,nsppol,spin
 81  character(len=fnlen) ::  seed_name  !seed names of files containing cg info used in case of MPI
 82  type(MPI_type),intent(in) :: mpi_enreg
 83 !arrays
 84  integer,intent(in) :: g1(3,nkpt,nntot),kg(3,mpw*mkmem),ngfft(18),npwarr(nkpt)
 85  integer,intent(in) :: iwav(mband,nkpt,nsppol)
 86  integer,intent(in) :: ovikp(nkpt,nntot)
 87  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
 88  real(dp),intent(out) :: cm1(2,mband,mband,nntot,nkpt,nsppol)
 89 
 90 !Local variables-------------------------------
 91 !scalars
 92  integer :: iband1,iband2,ierr,ig,ig1,ig1b,ig2,ig2b
 93  integer :: ig3,ig3b,igk1,igk2,igks1,igks2,ii,ikg,ikpt,ikpt1,ikpt2,imntot,index,intot,ios,ipw
 94  integer :: ispinor,isppol,iunit,me,n1,n2,n3,npoint,npoint2,npw_k,npw_k2
 95  integer :: nprocs,spaceComm
 96  integer,allocatable::indpwk(:,:),kg_k(:,:)
 97  integer,allocatable :: invpwk(:,:)   
 98  character(len=500) :: message
 99  character(len=fnlen) ::  cg_file  !file containing cg info used in case of MPI
100  logical::lfile
101  real(dp),allocatable :: cg_read(:,:) !to be used in case of MPI
102 
103 !************************************************************************
104 
105  write(message, '(a,a)' ) ch10,&
106 & '** mlwfovlp_pw : compute pw part of overlap'
107  call wrtout(std_out,  message,'COLL')
108  
109 !initialize flags
110  lfile=.false.
111 !mpi initialization
112  spaceComm=MPI_enreg%comm_cell
113  nprocs=xmpi_comm_size(spaceComm)
114  me=MPI_enreg%me_kpt
115 
116  if(nprocs>1) then
117    ABI_ALLOCATE(cg_read,(2,nspinor*mpw*mband))
118  end if
119 
120 
121 !****************compute intermediate quantities  (index, shifts) ******
122 !------------compute index for g points--------------------------------
123 !ig is a plane waves which belongs to the sphere ecut for ikpt (they
124 !are npwarr(ikpt))
125 !npoint is the position in the grid of planes waves
126 !(they are nfft)
127 !indpwk is a application ig-> npoint
128 !invpwk is not an application (some npoint have no ig corresponding)
129 !cg are ordered with npw_k !
130 !----------------------------------------------------------------------
131 !------------compute index for g points--------------------------------
132 !----------------------------------------------------------------------
133  write(message, '(a,a)' ) ch10,&
134 & '   first compute index for g-points'
135  call wrtout(std_out,  message,'COLL')
136 !
137 !Allocations
138  ABI_ALLOCATE(kg_k,(3,mpw))
139  ABI_ALLOCATE(indpwk,(nkpt,mpw))
140  ABI_ALLOCATE(invpwk,(nkpt,nfft))
141 !
142  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
143  invpwk=0
144  indpwk=0
145  kg_k=0
146  do isppol=1,1  !invpwk is not spin dependent
147 !  so we just do it once
148    ikg=0
149    do ikpt=1,nkpt
150 !    
151 !    MPI:cycle over k-points not treated by this node
152 !    
153      if ( ABS(MPI_enreg%proc_distrb(ikpt,1,isppol)-me)  /=0) CYCLE
154 
155 !    
156 !    write(std_out,*)'me',me,'ikpt',ikpt,'isppol',isppol
157      do npoint=1,nfft
158        if(invpwk(ikpt,npoint)/=0 )then
159          write(std_out,*) "error0 , invpwk is overwritten"
160          write(std_out,*) ikpt,npoint
161          MSG_ERROR("Aborting now")
162        end if
163      end do
164      npw_k=npwarr(ikpt)
165 !    write(std_out,*) ikpt,npw_k,nfft
166      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
167      do ig=1,npw_k
168        if(ig.gt.mpw) then
169          write(std_out,*)"error ig",ig,"greater than mpw ",mpw
170          MSG_ERROR("Aborting now")
171        end if
172        if(indpwk(ikpt,ig)/=0) then
173          write(std_out,*) "error, indpwk is overwritten"
174          write(std_out,*) ikpt,ig,indpwk(ikpt,ig)
175          MSG_ERROR("Aborting now")
176        end if
177        ig1=modulo(kg_k(1,ig),n1)
178        ig2=modulo(kg_k(2,ig),n2)
179        ig3=modulo(kg_k(3,ig),n3)
180        indpwk(ikpt,ig)=ig1+1+n1*(ig2+n2*ig3)
181        npoint=indpwk(ikpt,ig)
182        if(npoint.gt.nfft) then
183          MSG_ERROR("error npoint")
184        end if
185 !      write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint)
186        if(invpwk(ikpt,npoint)/=0) then
187          write(std_out,*) "error, invpwk is overwritten"
188          write(std_out,*) ikpt,ig,npoint,invpwk(ikpt,npoint)
189          MSG_ERROR("Aborting now")
190        end if
191        invpwk(ikpt,npoint)=ig
192 !      write(std_out,*)'ikpt,npoint,invpwk',ikpt,npoint,invpwk(ikpt,npoint)
193 !      if(ikpt.eq.1) write(std_out,*) "ig npoint",ig, npoint
194 !      write(std_out,*) "ikpt ig npoint",ikpt,ig, npoint
195      end do
196      ikg=ikg+npw_k
197 
198    end do !ikpt
199  end do !isppol
200 !write(std_out,*) "index for g points has been computed"
201 
202  call xmpi_barrier(spaceComm)
203  call xmpi_sum(invpwk,spaceComm,ierr)
204 
205 !----------------------------------------------------------------------
206 !------------test invpwk-----------------------------------------------
207 !----------------------------------------------------------------------
208 !write(std_out,*) "TEST INVPWK"
209 !ikpt=3
210 !isppol=1
211 !do ig=1,npwarr(ikpt)
212 !npoint=indpwk(ikpt,ig)
213 !write(std_out,*) "ig npoint    ",ig, npoint
214 !write(std_out,*) "ig npoint inv",invpwk(ikpt,npoint),npoint
215 !end do
216 !do ig3=1,n3
217 !do ig2=1,n2
218 !do ig1=1,n1
219 !npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1
220 !ig=invpwk(ikpt,npoint)
221 !!   if(ig/=0)  write(std_out,*) "ig npoint",ig, npoint
222 !end do
223 !end do
224 !end do
225 
226 
227 
228  
229 !
230 !Deallocate unused variables
231 !
232  ABI_DEALLOCATE(kg_k)
233  ABI_DEALLOCATE(indpwk)
234 
235 
236 !***********************************************************************
237 !**calculate overlap M_{mn}(k,b)=<\Psi_{k,m}|e^{-ibr}|\Psi_{k+b,n}>*****
238 !***********************************************************************
239  write(message, '(a,a)' ) ch10,&
240 & '   mlwfovlp_pw : compute overlaps '
241  call wrtout(std_out,  message,'COLL')
242  write(message, '(a,a)' ) ch10,&
243 & "     nkpt  nntot  mband "
244  call wrtout(std_out,  message,'COLL')
245  write(message, '(i6,2x,i6,2x,i6,2x,i6)' ) &
246 & nkpt,nntot,mband
247  call wrtout(std_out,  message,'COLL')
248  cm1=zero
249  write(message, '(a)' )  '  '
250  call wrtout(std_out,  message,'COLL')
251  do isppol=1,nsppol
252    if(spin.ne.0 .and. spin.ne.isppol) cycle
253    imntot=0
254    do ikpt1=1,nkpt
255 !    
256 !    MPI:cycle over k-points not treated by this node
257 !    
258      if ( ABS(MPI_enreg%proc_distrb(ikpt1,1,isppol)-me)  /=0) CYCLE
259 !    
260      write(message, '(a,i6,a,i6,a,i6)' ) &
261 &     '     Processor',me,' computes k-point',ikpt1,' and spin=',isppol
262      call wrtout(std_out,  message,'COLL')
263 !    write(std_out,*)trim(message)
264 
265      do intot=1,nntot
266        lfile=.false. !flag to know if this kpt will be read from a file, see below
267        imntot=imntot+1
268        ikpt2= ovikp(ikpt1,intot)
269 !      write(std_out,*)'me',me,'ikpt1',ikpt1,'ikpt2',ikpt2,'intot',intot,'isppol',isppol
270 
271 !      
272 !      MPI: if ikpt2 not found in this processor then
273 !      read info from an unformatted file
274 !      
275        if ( ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-me)  /=0) then
276          lfile=.true.
277          write(cg_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol 
278          iunit=1000+ikpt2+ikpt2*(isppol-1)
279          npw_k2=npwarr(ikpt2)
280          open (unit=iunit, file=cg_file,form='unformatted',status='old',iostat=ios)
281          if(ios /= 0) then
282            write(message,*) " mlwfovlp_pw: file",trim(cg_file), "not found"
283            MSG_ERROR(message)
284          end if
285 !        
286          do iband2=1,mband
287            do ipw=1,npw_k2*nspinor
288              index=ipw+(iband2-1)*npw_k2*nspinor
289              read(iunit) (cg_read(ii,index),ii=1,2)
290 !            if(me==0 .and. ikpt2==4)write(300,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
291 !            if(me==1 .and. ikpt2==4)write(301,*)'ipw,iband2,index',ipw,iband2,index,cg_read(:,index)
292            end do
293          end do
294          close(iunit)
295        end if
296 !      
297        npw_k=npwarr(ikpt1)
298        npw_k2=npwarr(ikpt2)
299        do ig3=1,n3
300          do ig2=1,n2
301            do ig1=1,n1
302 !            write(std_out,*) isppol,ikpt1,iband1,iband2,intot
303              npoint=ig1+(ig2-1)*n1+(ig3-1)*n2*n1
304              if(npoint.gt.nfft) then
305                write(std_out,*) "error npoint  b"
306                MSG_ERROR("Aborting now")
307              end if
308              ig1b=ig1+g1(1,ikpt1,intot)
309              ig2b=ig2+g1(2,ikpt1,intot)
310              ig3b=ig3+g1(3,ikpt1,intot)
311 !            write(std_out,*) ig1,ig2,ig3
312 !            write(std_out,*) ig1b,ig2b,ig3b
313              if(ig1b.lt.1) ig1b=ig1b+n1
314              if(ig2b.lt.1) ig2b=ig2b+n2
315              if(ig3b.lt.1) ig3b=ig3b+n3
316              if(ig1b.gt.n1) ig1b=ig1b-n1
317              if(ig2b.gt.n2) ig2b=ig2b-n2
318              if(ig3b.gt.n3) ig3b=ig3b-n3
319              npoint2=ig1b+(ig2b-1)*n1+(ig3b-1)*n2*n1
320              if(npoint2.gt.nfft) then
321                write(std_out,*)"error npoint  c"
322                MSG_ERROR("Aborting now")
323              end if
324              igk1=invpwk(ikpt1,npoint)
325              igk2=invpwk(ikpt2,npoint2) 
326              
327 !            if(intot==10) write(std_out,*)'Before igk1 and igk2',ikpt1,ikpt2,isppol
328 
329              if(igk1/=0.and.igk2/=0) then
330                do iband2=1,mband
331                  do iband1=1,mband
332                    do ispinor=1,nspinor
333 !                    igks1= (igk1*nspinor)-(nspinor-ispinor)
334 !                    igks2= (igk2*nspinor)-(nspinor-ispinor)
335                      igks1= igk1+ (ispinor-1)*npw_k
336                      igks2= igk2+ (ispinor-1)*npw_k2
337 
338 !                    Here the igks is to include the spinor component missing in igk
339                      if(lfile) index=igks2+npw_k2*nspinor*(iband2-1) !In case of MPI, see below
340 !                    
341 !                    If MPI sometimes the info was read from an unformatted file
342 !                    If that is the case lfile==.true.
343 !                    
344                      if(lfile) then
345                        cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ &
346 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index)&
347 &                       + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index)
348                        cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ &
349 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg_read(2,index)&
350 &                       - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg_read(1,index)
351 !                      
352                      else
353                        cm1(1,iband1,iband2,intot,ikpt1,isppol)=cm1(1,iband1,iband2,intot,ikpt1,isppol)+ &
354 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol))&
355 &                       + cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol))
356                        cm1(2,iband1,iband2,intot,ikpt1,isppol)=cm1(2,iband1,iband2,intot,ikpt1,isppol)+ &
357 &                       cg(1,igks1+iwav(iband1,ikpt1,isppol))*cg(2,igks2+iwav(iband2,ikpt2,isppol))&
358 &                       - cg(2,igks1+iwav(iband1,ikpt1,isppol))*cg(1,igks2+iwav(iband2,ikpt2,isppol))
359                      end if
360                    end do !ispinor
361                  end do ! iband1
362                end do ! iband2
363              end if
364            end do
365          end do
366        end do
367      end do ! intot
368    end do ! ikpt1
369  end do ! isppol
370 !
371 !Deallocations
372 !
373  ABI_DEALLOCATE(invpwk)
374  if(nprocs>1)  then
375    ABI_DEALLOCATE(cg_read)
376  end if
377 
378  end subroutine mlwfovlp_pw