TABLE OF CONTENTS


ABINIT/m_symkpt [ Modules ]

[ Top ] [ Modules ]

NAME

  m_symkpt

COPYRIGHT

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

TODO

  Move it to m_kpts

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_symkpt
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_sort
28  use m_krank
29  use m_numeric_tools
30  use m_time
31 
32  use m_fstrings,   only : sjoin, itoa
33 
34  implicit none
35 
36  private

ABINIT/mapkptsets [ Functions ]

[ Top ] [ Functions ]

NAME

 mapkptsets

FUNCTION

 given 2 input sets of kpts (1 and 2) find the symmetry operations that yield points in list 2 from a minimal set from list 1
 typical usage is to find k in a list from disk, to initialize the wfk in memory
 kin can be overcomplete etc... we just need to find _a_ solution
 bz2kin_smap follows the listkk convention with ik, isym, g0(1:3), itimrev

INPUTS

OUTPUT

 nkirred = number of irreducible k needed from list 1 (in)
 bz2ibz_smap = mapping of indices in list 2 with their irreducible origin in list 1,
   symop and timrev needed to transform them

SOURCE

670 subroutine mapkptsets(chksymbreak,gmet,k_in,nk_in,&
671 &   kbz,nkbz,nkirred,nsym,symrec,timrev,bz2kin_smap, comm)
672 
673 !Arguments -------------------------------
674 !scalars
675  integer,intent(in) :: chksymbreak,nkbz,nsym,timrev,comm
676  integer,intent(in) :: nk_in
677  integer,intent(out) :: nkirred
678 !arrays
679  integer,intent(in) :: symrec(3,3,nsym)
680  real(dp),intent(in) :: gmet(3,3),kbz(3,nkbz)
681  real(dp),intent(in) :: k_in(3,nk_in)
682  integer,intent(out) :: bz2kin_smap(nkbz,6)
683 
684 !Local variables -------------------------
685 !scalars
686  type(krank_t) :: krank
687  integer :: identi,ii,ikpt,ik_in,ikpt_found
688  integer :: isym,itim,jj,tident
689  !real(dp) :: cpu, gflops, wall
690  character(len=500) :: message
691 !arrays
692  real(dp) :: ksym(3),kpt1(3)
693 
694 ! *********************************************************************
695 
696  ABI_UNUSED(comm)
697  ABI_UNUSED(gmet)
698 
699  if (timrev/=1 .and. timrev/=0) then
700    write(message,'(a,i0)')' timrev should be 0 or 1, while it is equal to ',timrev
701    ABI_BUG(message)
702  end if
703 
704  ! Find the identity symmetry operation
705  identi = 1
706  tident = -1
707  if (nsym/=1) then
708    do isym=1,nsym
709      tident=1
710      do jj=1,3
711        if(symrec(jj,jj,isym)/=1)tident=0
712        do ii=1,3
713          if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
714        end do
715      end do
716      if(tident==1)then
717        identi=isym
718        exit
719      end if
720    end do
721    ABI_CHECK(tident == 1, 'Did not find the identity operation')
722  end if
723 
724  ! Initialize
725  bz2kin_smap = 0
726  do ikpt=1,nkbz
727    bz2kin_smap(ikpt, 1) = ikpt
728    bz2kin_smap(ikpt, 2) = 1
729    bz2kin_smap(ikpt, 3) = 1 ! We will use this as wtk_folded
730  end do
731 
732  ! Start krank
733  krank = krank_new(nkbz, kbz)
734 
735  ! Here begins the serious business
736  !call cwtime(cpu, wall, gflops, "start")
737 
738  ! If there is some possibility for a change
739  if(nkbz/=1 .and. (nsym/=1 .or. timrev==1) )then
740 
741    ! Examine whether the k point grid is symmetric or not
742    ! This check scales badly with nkbz hence it's disabled for dense meshes.
743    if (chksymbreak == 1 .and. nkbz < 40**3) then
744      do ikpt=1,nkbz
745        kpt1 = kbz(:,ikpt)
746 
747        do isym=1,nsym
748          do itim=0,timrev
749            ! Skip identity symmetry
750            if (isym==identi .and. itim==0) cycle
751 
752            ! Get the symmetric of the vector
753            do ii=1,3
754              ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
755                                    kpt1(2)*symrec(ii,2,isym)+&
756                                    kpt1(3)*symrec(ii,3,isym) )
757            end do
758 
759            !find this point
760            ikpt_found = krank%get_index(ksym)
761            !if (sum(abs(mod(ksym-kbz(:,ikpt_found),one)))>tol8) then
762            !  ABI_ERROR('Wrong k-point mapping found by krank')
763            !end if
764            !if k-point not found
765            if (ikpt_found < 0) then
766              write(message,'(3a,i4,2a,9i3,2a,i6,1a,3es16.6,6a)' )&
767              'Chksymbreak=1. It has been observed that the k point grid is not symmetric:',ch10,&
768              'for the symmetry number: ',isym,ch10,&
769              'with symrec= ',symrec(1:3,1:3,isym),ch10,&
770              'the symmetric of the k point number: ',ikpt,' with components: ',kpt1(:),ch10,&
771              'does not belong to the k point grid.',ch10,&
772              'Read the description of the input variable chksymbreak,',ch10,&
773              'You might switch it to zero, or change your k point grid to one that is symmetric.'
774              ABI_ERROR(message)
775            end if
776          end do ! itim
777        end do ! isym
778      end do ! ikpt
779    end if
780  end if ! End check on possibility of change
781  !call cwtime_report(" ibz", cpu, wall, gflops)
782 
783  ! Initialize
784  bz2kin_smap = 0
785 
786  ! HM: Here I invert the itim and isym loop to generate the same mapping as listkk
787  do itim=0,timrev
788    do isym=1,nsym
789 
790 !TODO: verify this inefficient use: sweep over isym=identity first to check which
791 !  k_in are actually directly present in kbz. Will not minimize the number of k_in we use, on the contrary
792 ! Now I loop over the points in the in list to find the mapping to the BZ
793    do ik_in=1,nk_in
794      kpt1 = k_in(:,ik_in)
795 
796        ! Get the symmetric of the vector
797        do ii=1,3
798          ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
799                                kpt1(2)*symrec(ii,2,isym)+&
800                                kpt1(3)*symrec(ii,3,isym) )
801        end do
802 
803        !find this point in the main set 2
804        ikpt_found = krank%get_index(ksym)
805 
806        if (ikpt_found < 0) cycle
807        ! if we already filled it, ignore new symmetric pre-image
808        if (bz2kin_smap(ikpt_found, 1) /= 0) cycle
809 
810        bz2kin_smap(ikpt_found,   1) = ik_in
811        bz2kin_smap(ikpt_found,   2) = isym
812        bz2kin_smap(ikpt_found, 3:5) = nint(kbz(:,ikpt_found)-ksym)
813        bz2kin_smap(ikpt_found,   6) = itim
814      end do
815 
816    end do
817  end do
818  !call cwtime_report(" map", cpu, wall, gflops)
819 
820  nkirred = 0
821  do ik_in=1,nk_in
822    ! did I end up using this ik_in?
823    if (any(bz2kin_smap(:,1) == ik_in)) then
824      nkirred = nkirred + 1
825    end if
826  end do
827 
828 ! check for redundant k in the kbz set
829  do ikpt=1, nkbz
830    if (bz2kin_smap(ikpt,1) /= 0) cycle
831 ! find index according to krank
832    ikpt_found = krank%get_index(kbz(:,ikpt))
833 ! if I am not my own hash image, associate the smap data from my image
834    if (ikpt_found /= ikpt) then
835      bz2kin_smap(ikpt,:) = bz2kin_smap(ikpt_found,:)
836    end if
837  end do
838 
839  call krank%free()
840 
841  !Here I make a check if the mapping was sucessfull
842  !might exit less brutally and allow for error catching by the caller...
843  if (any(bz2kin_smap(:,1) == 0)) then
844 !print *, 'nkirred ', nkirred
845 !do ikpt_found=1, nkbz
846 !print *, bz2kin_smap(ikpt_found,:), kbz(:,ikpt_found)
847 !end do
848 !print *, 'k_in = '
849 !do ik_in=1, nk_in
850 !print *, k_in(:,ik_in)
851 !end do
852    ABI_ERROR('Could not find mapping k-point sets')
853  end if
854 
855  !do ikpt=1,nkbz
856  !  write(*,*) ikpt, ibz2bz(ikpt), bz2kin_smap(ikpt,1)
857  !end do
858 
859 end subroutine mapkptsets

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).

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.
 comm=MPI communicator.

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
 bz2ibz_smap(nkbz, 6)= Mapping BZ --> IBZ.

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

TODO

  Bad scaling wrt nkbz. Should try to MPI parallelize or implement more efficient algorithm

SOURCE

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

ABINIT/symkpt_new [ Functions ]

[ Top ] [ Functions ]

NAME

 symkpt_new

FUNCTION

 Same routine as above but with an algorithm with better scalling than before.
 From a few tests it produces the same IBZ as before but avoids computing the lengths and sorting.
 Instead it uses the krank datatype to map k-points onto each other.

INPUTS

OUTPUT

SOURCE

424 subroutine symkpt_new(chksymbreak,gmet,ibz2bz,iout,kbz,nkbz,nkibz,nsym,symrec,timrev,bz2ibz_smap, comm)
425 
426 !Arguments -------------------------------
427 !scalars
428  integer,intent(in) :: chksymbreak,iout,nkbz,nsym,timrev,comm
429  integer,intent(out) :: nkibz
430 !arrays
431  integer,intent(in) :: symrec(3,3,nsym)
432  integer,intent(out) :: ibz2bz(nkbz)
433  real(dp),intent(in) :: gmet(3,3),kbz(3,nkbz)
434  integer,intent(out) :: bz2ibz_smap(6,nkbz)
435 
436 !Local variables -------------------------
437 !scalars
438  type(krank_t) :: krank
439  integer :: identi,ii,ikpt,ikibz,ikpt_found
440  integer :: isym,itim,jj,nkpout,tident
441  !real(dp) :: cpu, gflops, wall
442  character(len=500) :: msg
443 !arrays
444  real(dp) :: ksym(3),kpt1(3)
445 
446 ! *********************************************************************
447 
448  ABI_UNUSED(comm)
449  ABI_UNUSED(gmet)
450 
451  if (timrev/=1 .and. timrev/=0) then
452    write(msg,'(a,i0)')' timrev should be 0 or 1, while it is equal to ',timrev
453    ABI_BUG(msg)
454  end if
455 
456  ! Find the identity symmetry operation
457  identi = 1
458  tident = -1
459  if (nsym/=1) then
460    do isym=1,nsym
461      tident=1
462      do jj=1,3
463        if(symrec(jj,jj,isym)/=1)tident=0
464        do ii=1,3
465          if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
466        end do
467      end do
468      if(tident==1)then
469        identi=isym
470        !call wrtout(std_out, sjoin(' symkpt_new: found identity, with number: ',itoa(identi)))
471        exit
472      end if
473    end do
474    ABI_CHECK(tident == 1, 'Did not find the identity operation')
475  end if
476 
477  ! Initialize
478  ibz2bz = 0
479  bz2ibz_smap = 0
480  do ikpt=1,nkbz
481    bz2ibz_smap(1, ikpt) = ikpt
482    bz2ibz_smap(2, ikpt) = 1
483    bz2ibz_smap(4, ikpt) = 1 ! We will use this as wtk_folded
484  end do
485 
486  ! Start krank
487  krank = krank_new(nkbz, kbz)
488 
489  ! Here begins the serious business
490  !call cwtime(cpu, wall, gflops, "start")
491 
492  ! If there is some possibility for a change
493  if(nkbz/=1 .and. (nsym/=1 .or. timrev==1) )then
494 
495    ! Examine whether the k point grid is symmetric or not
496    ! This check scales badly with nkbz hence it's disabled for dense meshes.
497    if (chksymbreak == 1 .and. nkbz < 40**3) then
498      do ikpt=1,nkbz
499        kpt1 = kbz(:,ikpt)
500 
501        do isym=1,nsym
502          do itim=0,timrev
503            ! Skip identity symmetry
504            if (isym==identi .and. itim==0) cycle
505 
506            ! Get the symmetric of the vector
507            do ii=1,3
508              ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
509                                    kpt1(2)*symrec(ii,2,isym)+&
510                                    kpt1(3)*symrec(ii,3,isym) )
511            end do
512 
513            !find this point
514            ikpt_found = krank%get_index(ksym)
515            !if (sum(abs(mod(ksym-kbz(:,ikpt_found),one)))>tol8) then
516            !  ABI_ERROR('Wrong k-point mapping found by krank')
517            !end if
518            !if k-point not found
519            if (ikpt_found < 0) then
520              write(msg,'(3a,i4,2a,9i3,2a,i6,1a,3es16.6,6a)' )&
521              'Chksymbreak=1. It has been observed that the k point grid is not symmetric:',ch10,&
522              'for the symmetry number: ',isym,ch10,&
523              'with symrec= ',symrec(1:3,1:3,isym),ch10,&
524              'the symmetric of the k point number: ',ikpt,' with components: ',kpt1(:),ch10,&
525              'does not belong to the k point grid.',ch10,&
526              'Read the description of the input variable chksymbreak,',ch10,&
527              'You might switch it to zero, or change your k point grid to one that is symmetric.'
528              ABI_ERROR(msg)
529            end if
530          end do ! itim
531        end do ! isym
532      end do ! ikpt
533    end if
534 
535    ! Here I generate the IBZ
536    ikpt_loop: do ikpt=1,nkbz
537 
538      if (bz2ibz_smap(4, ikpt)==0) cycle
539      kpt1 = kbz(:,ikpt)
540 
541      ! MG Dec 16 2018, Invert isym, itim loop to be consistent with listkk and GW routines
542      ! Should always use this convention when applying symmetry operations in k-space.
543      ! TODO: Postponed to v9 because it won't be possible to read old WFK files.
544      do isym=1,nsym
545        do itim=0,timrev
546          ! Skip identity symmetry
547          if (isym==identi .and. itim==0) cycle
548 
549          ! Get the symmetric of the vector
550          do ii=1,3
551            ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
552                                  kpt1(2)*symrec(ii,2,isym)+&
553                                  kpt1(3)*symrec(ii,3,isym) )
554          end do
555 
556          !find this point
557          ikpt_found = krank%get_index(ksym)
558          !if k-point not found just cycle
559          if (ikpt_found < 0) cycle
560          !if (sum(abs(mod(ksym-kbz(:,ikpt_found),one)))>tol8) then
561          !  ABI_ERROR('Wrong k-point mapping found by krank')
562          !end if
563          if (ikpt_found >= ikpt) cycle
564          bz2ibz_smap(:3, ikpt)  = [ikpt_found, isym, itim]
565          bz2ibz_smap(4,ikpt) = bz2ibz_smap(4,ikpt) + bz2ibz_smap(4,ikpt_found)
566          bz2ibz_smap(4,ikpt_found) = 0
567        end do ! itim
568      end do ! isym
569    end do ikpt_loop! ikpt
570 
571  end if ! End check on possibility of change
572  !call cwtime_report(" ibz", cpu, wall, gflops)
573 
574  nkibz = 0
575  do ikpt=1,nkbz
576    ikibz = bz2ibz_smap(1,ikpt)
577    if (ikibz /= ikpt) cycle
578    nkibz = nkibz + 1
579    ibz2bz(nkibz) = ikpt
580  end do
581 
582  !do ikpt=1,nkbz
583  !  write(*,*) ikpt, ibz2bz(ikpt), bz2ibz_smap(1,ikpt)
584  !end do
585 
586  ! Initialize again
587  bz2ibz_smap = 0
588 
589  ! Now I loop again over the points in the IBZ to find the mapping to the BZ
590  do ikpt=1,nkibz
591    ikibz = ibz2bz(ikpt)
592    kpt1 = kbz(:,ikibz)
593 
594    ! HM: Here I invert the itim and isym loop to generate the same mapping as listkk
595    do itim=0,timrev
596      do isym=1,nsym
597        ! Get the symmetric of the vector
598        do ii=1,3
599          ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
600                                kpt1(2)*symrec(ii,2,isym)+&
601                                kpt1(3)*symrec(ii,3,isym) )
602        end do
603 
604        ! Find this point
605        ikpt_found = krank%get_index(ksym)
606        if (ikpt_found < 0) cycle
607        if (bz2ibz_smap(1, ikpt_found) /= 0) cycle
608        bz2ibz_smap(:3, ikpt_found) = [ikpt, isym, itim]
609        bz2ibz_smap(4:, ikpt_found) = nint(kbz(:,ikpt_found)-ksym)
610 
611        ! TODO: Use same conventions as in listkk but must propagate the changes!
612        !bz2ibz_smap(1:2, ikpt_found) = [ikpt, isym]
613        !bz2ibz_smap(3:5, ikpt_found) = nint(-ksym + kbz(:,ikpt_found)  + tol12)
614        !bz2ibz_smap(6, ikpt_found) = itim
615      end do
616    end do
617  end do
618  !call cwtime_report(" map", cpu, wall, gflops)
619 
620  call krank%free()
621 
622  ! Check whether the mapping was sucessfull
623  if (any(bz2ibz_smap(1, :) == 0)) then
624    ABI_ERROR('Could not initial mapping BZ to IBZ. Perhaps your grid breaks the symmetry of the lattice!')
625  end if
626 
627  !do ikpt=1,nkbz
628  !  write(*,*) ikpt, ibz2bz(ikpt), bz2ibz_smap(1,ikpt)
629  !end do
630 
631  if(iout/=0)then
632    if(nkbz/=nkibz)then
633      write(msg, '(a,a,a,i6,a)' )&
634      ' symkpt_new : the number of k-points, thanks to the symmetries,',ch10,' is reduced to',nkibz,' .'
635      call wrtout(iout,msg)
636      if(iout/=std_out) call wrtout(std_out,msg)
637 
638      nkpout=nkibz
639    else
640      write(msg, '(a)' )' symkpt_new : not enough symmetry to change the number of k points.'
641      call wrtout(iout,msg)
642      if (iout/=std_out) call wrtout(std_out,msg)
643    end if
644  end if
645 
646 end subroutine symkpt_new