TABLE OF CONTENTS


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

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