TABLE OF CONTENTS


ABINIT/m_symkpt [ Modules ]

[ Top ] [ Modules ]

NAME

  m_symkpt

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

TODO

  Move it to m_kpts

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 module m_symkpt
31 
32  implicit none
33 
34  private

ABINIT/symkpt [ Functions ]

[ Top ] [ Functions ]

NAME

 symkpt

FUNCTION

 Determines the weights of the k-points for sampling the Brillouin Zone, starting from a first set
 of weights wtk, and folding it to a new set, by taking into account the symmetries described
 by symrec, and eventually the time-reversal symmetry.
 Also compute the number of k points in the reduced set
 This routine is also used for sampling the q vectors in the Brillouin zone for the computation
 of thermodynamical properties (from the routine thm9).

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (XG,LSI)
 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

 chksymbreak= if 1, will check whether the k point grid is symmetric, and stop if not.
 gmet(3,3)=reciprocal space metric (bohr**-2).
 iout=if non-zero, output the new number of kpoints on unit iout
 kbz(3,nkbz)= k vectors in the BZ.
 nkbz = number of k-points whose weights are wtk
 nsym=number of space group symmetries
 symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
 timrev: if 1, the time reversal operation has to be taken into account
         if 0, no time reversal symmetry.
 wtk(nkbz)=weight assigned to each k point.

OUTPUT

 ibz2bz(nkbz)=non-symmetrized indices of the k-points a.k.a. ibz2bz mapping
   The correspondence beween the iq_ibz point in IBZ and the iq_bz point in the full BZ is obtained via:

       do ik_ibz=1,nkibz
         ik_bz = ibz2bz(ik_ibz)
       end do

 nkibz = number of k-points in the irreducible set
 wtk_folded(nkbz)=weight assigned to each k point, taking into account the symmetries

NOTES

 The decomposition of the symmetry group in its primitives might speed up the execution.
 The output variables are stored only in the range 1:nkbz

PARENTS

      dfpt_looppert,elphon,ep_setupqpt,get_npert_rbz,getkgrid,harmonic_thermo
      m_bz_mesh,m_ifc,m_sigmaph

CHILDREN

      sort_dp,wrtout

SOURCE

 98 subroutine symkpt(chksymbreak,gmet,ibz2bz,iout,kbz,nkbz,nkibz,nsym,&
 99 & symrec,timrev,wtk,wtk_folded)
100 
101  use defs_basis
102  use m_abicore
103  use m_errors
104  use m_sort
105 
106 !This section has been created automatically by the script Abilint (TD).
107 !Do not modify the following lines by hand.
108 #undef ABI_FUNC
109 #define ABI_FUNC 'symkpt'
110 !End of the abilint section
111 
112  implicit none
113 
114 !Arguments -------------------------------
115 !scalars
116  integer,intent(in) :: chksymbreak,iout,nkbz,nsym,timrev
117  integer,intent(out) :: nkibz
118 !arrays
119  integer,intent(in) :: symrec(3,3,nsym)
120  integer,intent(inout) :: ibz2bz(nkbz) !vz_i
121  real(dp),intent(in) :: gmet(3,3),kbz(3,nkbz),wtk(nkbz)
122  real(dp),intent(out) :: wtk_folded(nkbz)
123 
124 !Local variables -------------------------
125 !scalars
126  integer :: identi,ii,ikpt,ikpt2,ind_ikpt,ind_ikpt2
127  integer :: ikpt_current_length,isym,itim,jj,nkpout,quit,tident
128  real(dp) :: difk,difk1,difk2,difk3,length2trial,reduce,reduce1,reduce2,reduce3
129  character(len=500) :: message
130 !arrays
131  integer,allocatable :: list(:)
132  real(dp) :: gmetkpt(3),ksym(3)
133  real(dp),allocatable :: length2(:)
134 
135 ! *********************************************************************
136 
137  if (timrev/=1 .and. timrev/=0) then
138    write(message,'(a,i0)')' timrev should be 0 or 1, while it is equal to ',timrev
139    MSG_BUG(message)
140  end if
141 
142  if(nsym/=1)then
143 !  Find the identity symmetry operation
144    do isym=1,nsym
145      tident=1
146      do jj=1,3
147        if(symrec(jj,jj,isym)/=1)tident=0
148        do ii=1,3
149          if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
150        end do
151      end do
152      if(tident==1)then
153        identi=isym
154        !write(message, '(a,i0)' )' symkpt : found identity, with number',identi
155        !call wrtout(std_out,message,'COLL')
156        exit
157      end if
158    end do
159    if(tident==0)then
160      MSG_BUG(' Did not find the identity operation')
161    end if
162  else
163    identi=1
164  end if
165 
166 !Initialise the wtk_folded array using the wtk array :
167  do ikpt=1,nkbz
168    wtk_folded(ikpt)=wtk(ikpt)
169  end do
170 
171 !Here begins the serious business
172 
173 !If there is some possibility for a change (otherwise, wtk_folded is
174 !correctly initialized to give no change) :
175  if(nkbz/=1 .and. (nsym/=1 .or. timrev==1) )then
176 
177 !  Store the length of vectors, but take into account umklapp
178 !  processes by selecting the smallest length of all symmetric vectors
179    ABI_ALLOCATE(length2,(nkbz))
180 
181 
182 !MG FIXME:
183 ! Here there's a possible problem with the order of symmetries because
184 ! in listkk, time-reversal is the outermost loop. This can create inconsistencies in the symmetry tables.
185    do ikpt=1,nkbz
186      do isym=1,nsym
187        do itim=1,(1-2*timrev),-2
188 !        Get the symmetric of the vector
189          do ii=1,3
190            ksym(ii)=itim*( kbz(1,ikpt)*symrec(ii,1,isym)&
191 &           +kbz(2,ikpt)*symrec(ii,2,isym)&
192 &           +kbz(3,ikpt)*symrec(ii,3,isym) )
193            ksym(ii)=ksym(ii)-anint(ksym(ii)+tol8*half)
194          end do
195          gmetkpt(:)=gmet(:,1)*ksym(1)+gmet(:,2)*ksym(2)+gmet(:,3)*ksym(3)
196          length2trial=ksym(1)*gmetkpt(1)+ksym(2)*gmetkpt(2)+ksym(3)*gmetkpt(3)
197          if(isym==1 .and. itim==1)then
198            length2(ikpt)=length2trial
199          else
200            if(length2(ikpt)>length2trial)length2(ikpt)=length2trial
201          end if
202        end do
203      end do
204    end do
205 
206 !  Sort the lengths
207    ABI_ALLOCATE(list,(nkbz))
208    list(:)=(/ (ikpt,ikpt=1,nkbz) /)
209    call sort_dp(nkbz,length2,list,tol14)
210 !  do ikpt=1,nkbz
211 !  write(std_out,*)ikpt,length2(ikpt),list(ikpt),kbz(1:3,list(ikpt))
212 !  end do
213 
214 !  Examine whether the k point grid is symmetric or not
215    if(chksymbreak==1)then
216      ikpt_current_length=1
217 !    Loop on all k points
218      do ikpt=1,nkbz
219        ind_ikpt=list(ikpt)
220 !      Keep track of the current length, to avoid doing needless comparisons
221        if(length2(ikpt)-length2(ikpt_current_length)>tol8)then
222          ikpt_current_length=ikpt
223        end if
224 
225        do isym=1,nsym
226          do itim=1,(1-2*timrev),-2
227            if(isym/=identi .or. itim/=1 )then
228 
229 !            Get the symmetric of the vector
230              do ii=1,3
231                ksym(ii)=itim*( kbz(1,ind_ikpt)*symrec(ii,1,isym)&
232 &               +kbz(2,ind_ikpt)*symrec(ii,2,isym)&
233 &               +kbz(3,ind_ikpt)*symrec(ii,3,isym) )
234              end do
235 
236 !            Search over k-points with the same length, to find whether there is a connecting symmetry operation
237              quit=0
238              do ikpt2=ikpt_current_length,nkbz
239 !              The next line skip all ikpt2 vectors, as soon as one becomes larger than length2(ikpt)
240 !              Indeed, one is already supposed to have found a symmetric k point before this happens ...
241                if(length2(ikpt2)-length2(ikpt)>tol8)exit
242 !              Ordered index
243                ind_ikpt2=list(ikpt2)
244                difk1= ksym(1)-kbz(1,ind_ikpt2)
245                reduce1=difk1-anint(difk1)
246                difk2= ksym(2)-kbz(2,ind_ikpt2)
247                reduce2=difk2-anint(difk2)
248                difk3= ksym(3)-kbz(3,ind_ikpt2)
249                reduce3=difk3-anint(difk3)
250                if(abs(reduce1)+abs(reduce2)+abs(reduce3)<tol8)then
251 !                The symmetric was found
252                  quit=1
253                  exit
254                end if
255              end do
256              if(quit==0)then
257                write(message,'(3a,i4,2a,9i3,2a,i6,1a,3es16.6,6a)' )&
258 &               'Chksymbreak=1 . It has been observed that the k point grid is not symmetric :',ch10,&
259 &               'for the symmetry number ',isym,ch10,&
260 &               'with symrec=',symrec(1:3,1:3,isym),ch10,&
261 &               'the symmetric of the k point number ',ind_ikpt2,' with components', kbz(1:3,ind_ikpt2),ch10,&
262 &               'does not belong to the k point grid.',ch10,&
263 &               'Read the description of the input variable chksymbreak,',ch10,&
264 &               'You might switch it to zero, or change your k point grid to one that is symmetric.'
265                MSG_ERROR(message)
266              end if
267 
268            end if ! End condition of non-identity symmetry
269          end do ! itim
270        end do ! isym
271 
272      end do ! ikpt
273    end if ! chksymbreak==1
274 
275 
276 !  Eliminate the k points that are symmetric of another one
277    do ikpt=1,nkbz-1
278 
279 !    Ordered index
280      ind_ikpt=list(ikpt)
281 
282 !    Not worth to examine a k point that is a symmetric of another,
283 !    which is the case if its weight has been set to 0 by previous folding
284      if(wtk_folded(ind_ikpt)<tol16)cycle
285 
286 !    Loop on the remaining k-points
287      do ikpt2=ikpt+1,nkbz
288 
289 !      The next line eliminates pairs of vectors that differs by their length.
290 !      Moreover, since the list is ordered according to the length,
291 !      one can skip all other ikpt2 vectors, as soon as one becomes larger than length2(ikpt)
292        if(length2(ikpt2)-length2(ikpt)>tol8)exit
293 
294 !      Ordered index
295        ind_ikpt2=list(ikpt2)
296 
297 !      If the second vector is already empty, no interest to treat it
298        if(wtk_folded(ind_ikpt2)<tol16)cycle
299 
300        quit=0
301        do isym=1,nsym
302          do itim=1,(1-2*timrev),-2
303            if(isym/=identi .or. itim/=1 )then
304 
305 !            Get the symmetric of the vector
306              do ii=1,3
307                ksym(ii)=itim*( kbz(1,ind_ikpt)*symrec(ii,1,isym)&
308 &               +kbz(2,ind_ikpt)*symrec(ii,2,isym)&
309 &               +kbz(3,ind_ikpt)*symrec(ii,3,isym) )
310              end do
311 
312 !            The do-loop was expanded to speed up the execution
313              difk= ksym(1)-kbz(1,ind_ikpt2)
314              reduce=difk-anint(difk)
315              if(abs(reduce)>tol8)cycle
316              difk= ksym(2)-kbz(2,ind_ikpt2)
317              reduce=difk-anint(difk)
318              if(abs(reduce)>tol8)cycle
319              difk= ksym(3)-kbz(3,ind_ikpt2)
320              reduce=difk-anint(difk)
321              if(abs(reduce)>tol8)cycle
322 
323 !            Here, have successfully found a symmetrical k-vector
324 !            Assign all the weight of the k-vector to its symmetrical
325              wtk_folded(ind_ikpt)=wtk_folded(ind_ikpt)+wtk_folded(ind_ikpt2)
326              wtk_folded(ind_ikpt2)=0._dp
327 
328 !            Go to the next ikpt2 if the symmetric was found
329              quit=1
330              exit
331 
332            end if ! End condition of non-identity symmetry
333          end do ! End loop on itim
334 
335          if(quit==1)exit
336        end do !  End loop on isym
337      end do ! End secondary loop over k-points
338    end do ! End primary loop over k-points
339 
340    ABI_DEALLOCATE(length2)
341    ABI_DEALLOCATE(list)
342  end if ! End check on possibility of change
343 
344 !Create the indexing array ibz2bz
345  nkibz=0
346  do ikpt=1,nkbz
347    if(wtk_folded(ikpt)>tol8)then
348      nkibz=nkibz+1
349      ibz2bz(nkibz)=ikpt
350    end if
351  end do
352 
353  if(iout/=0)then
354    if(nkbz/=nkibz)then
355      write(message, '(a,a,a,i6,a)' )&
356 &     ' symkpt : the number of k-points, thanks to the symmetries,',ch10,&
357 &     ' is reduced to',nkibz,' .'
358      call wrtout(iout,message,'COLL')
359      if(iout/=std_out) then
360        call wrtout(std_out,message,'COLL')
361      end if
362 
363      nkpout=nkibz
364      !if(nkibz>80)then
365      !  write(message,'(a)' )' greater than 80, so only write 20 of them '
366      !  call wrtout(std_out,message,'COLL')
367      !  nkpout=20
368      !end if
369      !do ii=1,nkpout
370      !  write(message, '(1x,i2,a2,3es16.8)' ) ii,') ',kbz(1:3,ibz2bz(ii))
371      !  call wrtout(std_out,message,'COLL')
372      !end do
373 
374 !    DEBUG
375 !    write(message, '(a)' )'   Here are the new weights :'
376 !    call wrtout(std_out,message,'COLL')
377 !    do ikpt=1,nkbz,6
378 !    write(message, '(6f12.6)' ) wtk_folded(ikpt:min(nkbz,ikpt+5))
379 !    call wrtout(std_out,message,'COLL')
380 !    end do
381 !    ENDDEBUG
382    else
383      write(message, '(a)' )' symkpt : not enough symmetry to change the number of k points.'
384      call wrtout(iout,message,'COLL')
385      if(iout/=std_out) then
386        call wrtout(std_out,message,'COLL')
387      end if
388    end if
389  end if
390 
391 !DEBUG
392 !write(std_out,*)' exit symkpt '
393 !write(std_out,*)' nkbz,nsym',nkbz,nsym
394 !write(std_out,*)' wtk',wtk
395 !write(std_out,*)' timrev',timrev
396 !if(timrev==0)stop
397 !if(option==1)stop
398 !ENDDEBUG
399 
400 end subroutine symkpt