TABLE OF CONTENTS
ABINIT/m_symkpt [ Modules ]
NAME
m_symkpt
COPYRIGHT
Copyright (C) 1999-2024 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 ]
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 ]
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 ]
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