TABLE OF CONTENTS
ABINIT/listkk [ Functions ]
NAME
listkk
FUNCTION
Given a list of nkpt1 initial k points kptns1 and a list of nkpt2 final k points kptns2, associates each final k pt with a "closest" initial k point (or symmetric thereof, also taking possible umklapp) as determined by a metric gmet, that commutes with the symmetry operations. The algorithm does not scale as nkpt1 times nkpt2, thanks to the ordering of the kptns1 and kptns2 vectors according to their lengths, and comparison first between vectors of similar lengths. Returns indirect indexing list indkk.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR) 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
gmet(3,3)=reciprocal space metric (bohr^-2) kptns1(3,nkpt1)=list of initial k points (reduced coordinates) kptns2(3,nkpt2)=list of final k points nkpt1=number of initial k points nkpt2=number of final k points nsym=number of symmetry elements in space group sppoldbl=if 1, no spin-polarisation doubling if 2, spin-polarisation doubling using symafm symafm(nsym)=(anti)ferromagnetic part of symmetry operations symmat(3,3,nsym)=symmetry operations (symrel or symrec, depending on value of use_symrec timrev=1 if the use of time-reversal is allowed; 0 otherwise use_symrec: if present and true, symmat assumed to be symrec, otherwise assumed to be symrel (default)
OUTPUT
dksqmax=maximal value of the norm**2 of the difference between
a kpt2 vector and the closest k-point found from the kptns1 set, using symmetries.
indkk(nkpt2*sppoldbl,6)=describe k point number of kpt1 that allows to
generate wavefunctions closest to given kpt2
if sppoldbl=2, use symafm to generate spin down wfs from spin up wfs
indkk(:,1)=k point number of kptns1
indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
(if 0, means no symmetry operation, equivalent to identity )
indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
to give kpt1b, that is the closest to kpt2.
indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
NOTES
The tolerances tol12 and tol8 aims at giving a machine-independent ordering. (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f) The tolerance tol12 is used for each component of the k vectors, and for the length of the vectors while the tolerance tol8 is used for the comparison of the squared lengths of the separate vectors.
PARENTS
initberry,initorbmag,inwffil,m_dvdb,m_ebands,m_eprenorms,m_exc_diago m_fock,m_fstab,m_haydock,m_ifc,m_kpts,m_phgamma,m_sigmaph,mlwfovlp_qp
CHILDREN
sort_dp,timab
SOURCE
68 #if defined HAVE_CONFIG_H 69 #include "config.h" 70 #endif 71 72 #include "abi_common.h" 73 74 75 subroutine listkk(dksqmax,gmet,indkk,kptns1,kptns2,nkpt1,nkpt2,nsym,& 76 & sppoldbl,symafm,symmat,timrev,use_symrec) 77 78 use defs_basis 79 use m_errors 80 use m_profiling_abi 81 use m_sort 82 83 !This section has been created automatically by the script Abilint (TD). 84 !Do not modify the following lines by hand. 85 #undef ABI_FUNC 86 #define ABI_FUNC 'listkk' 87 use interfaces_18_timing 88 !End of the abilint section 89 90 implicit none 91 92 !Arguments ------------------------------------ 93 !scalars 94 integer,intent(in) :: nkpt1,nkpt2,nsym,sppoldbl,timrev 95 real(dp),intent(out) :: dksqmax 96 logical,optional,intent(in) :: use_symrec 97 !arrays 98 integer,intent(in) :: symafm(nsym),symmat(3,3,nsym) 99 integer,intent(out) :: indkk(nkpt2*sppoldbl,6) 100 real(dp),intent(in) :: gmet(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2) 101 102 !Local variables------------------------------- 103 !scalars 104 integer :: l3,ig1,ig2,ig3,ii,ikpg1,ikpt1,ikpt2,ikpt2_done 105 integer :: ilarger,ismaller,itrial 106 integer :: isppol,isym,itimrev,jkpt1,jsym,jtime,limit 107 integer :: nsym_used,timrev_used,usesym 108 real(dp) :: dksq,dksqmn,lk2,llarger,ldiff,lsmaller,ltrial,min_l 109 character(len=500) :: message 110 !arrays 111 integer :: dkint(3),jdkint(3),k1int(3),k2int(3) 112 integer, allocatable :: isort(:) 113 real(dp) :: tsec(2) 114 real(dp) :: dk(3),kpg1(3),kpt1a(3),k1(3),k2(3) 115 !real(dp) :: kasq,ka(3) 116 real(dp),allocatable :: lkpg1(:),lkpg1_sorted(:) 117 118 ! ************************************************************************* 119 120 !write(std_out,*)' listkk : nkpt1,nkpt2,nsym=',nkpt1,nkpt2,nsym 121 call timab(1021,1,tsec) 122 123 if(sppoldbl<1 .or. sppoldbl>2)then 124 write(message, '(a,i4,3a)' )& 125 & 'The value of sppoldbl is',sppoldbl,',',ch10,& 126 & 'but it should be either 1 or 2.' 127 MSG_BUG(message) 128 end if 129 130 !When usesym=0, the old way of converting the wavefunctions (without 131 !using the symmetries), is recovered. 132 usesym=1 133 134 nsym_used=nsym 135 timrev_used=timrev 136 if(usesym==0)nsym_used=1 137 if(usesym==0)timrev_used=0 138 139 !Precompute the length of the kpt1 vectors, also taking into account 140 !possible umpklapp vectors 141 limit=1 ; l3 = (2*limit+1)**3 142 ABI_ALLOCATE(lkpg1,(l3*nkpt1)) 143 ABI_ALLOCATE(lkpg1_sorted,(l3*nkpt1)) 144 ABI_ALLOCATE(isort,(l3*nkpt1)) 145 !write(std_out,*)' List of kpt1 vectors ' 146 !write(std_out,*)' Length of the kpt1 vectors :' 147 148 do ikpt1=1,nkpt1 149 k1(:)=kptns1(:,ikpt1) 150 ! write(std_out,*)ikpt1,k1(:) 151 k1int(:)=nint(k1(:)+tol12) 152 k1(:)=k1(:)-k1int(:) 153 do ig1=-limit,limit 154 kpg1(1)=k1(1)+ig1 155 do ig2=-limit,limit 156 kpg1(2)=k1(2)+ig2 157 do ig3=-limit,limit 158 kpg1(3)=k1(3)+ig3 159 160 ikpg1=ig1+limit+1 + (2*limit+1)*(ig2+limit) + (2*limit+1)**2*(ig3+limit) + l3*(ikpt1-1) 161 ! Compute the norm of the vector (also taking into account possible umklapp) 162 lkpg1(ikpg1)=sqrt(gmet(1,1)*kpg1(1)**2+gmet(2,2)*kpg1(2)**2+& 163 & gmet(3,3)*kpg1(3)**2+two*(gmet(2,1)*kpg1(2)*kpg1(1)+& 164 & gmet(3,2)*kpg1(3)*kpg1(2)+gmet(3,1)*kpg1(3)*kpg1(1))) 165 lkpg1_sorted(ikpg1)=lkpg1(ikpg1) 166 isort(ikpg1)=ikpg1 167 ! write(std_out,*)' ikpt1,ig1,ig2,ig3,lkpg1=',ikpt1,ig1,ig2,ig3,lkpg1(ikpg1) 168 end do 169 end do 170 end do 171 end do 172 173 call sort_dp( l3*nkpt1,lkpg1_sorted,isort,tol12) 174 175 !DEBUG 176 !write(std_out,*)' listkk : output list of kpt1 for checking purposes ' 177 !write(std_out,*)' ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii)) ' 178 !do ii=1,l3*nkpt1 179 !ikpt1=(isort(ii)-1)/l3+1 180 !write(std_out,*)ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii)) 181 !enddo 182 !stop 183 !ENDDEBUG 184 185 dksqmax=zero 186 do isppol=1,sppoldbl 187 do ikpt2=1,nkpt2 188 189 ikpt2_done=0 190 ! Precompute the length of the kpt2 vector, with the Umklapp vector such that it is the closest to the Gamma point 191 k2(:)=kptns2(:,ikpt2) 192 k2int(:)=nint(k2(:)+tol12) 193 k2(:)=k2(:)-k2int(:) 194 lk2=sqrt(gmet(1,1)*k2(1)**2+gmet(2,2)*k2(2)**2+& 195 & gmet(3,3)*k2(3)**2+two*(gmet(2,1)*k2(2)*k2(1)+& 196 & gmet(3,2)*k2(3)*k2(2)+gmet(3,1)*k2(3)*k2(1))) 197 198 ! DEBUG 199 ! write(std_out, '(a,i4,7es16.6)' )' listkk : ikpt2,kptns2(:,ikpt2),k2(:),lk2=',ikpt2,kptns2(:,ikpt2),k2(:),lk2 200 ! if(ikpt2/=17)cycle 201 ! ENDDEBUG 202 203 ! Find the kpt1 vector whose length is the most similar to the length of lk2 204 ! up to a tolerance. Use a bissection algorithm. 205 ismaller=0 ; lsmaller=zero 206 ilarger=l3*nkpt1+1 ; llarger=huge(one) 207 ! This loop should never reach l3*nkpt1, since this is a bissection algorithm 208 do ii=1,l3*nkpt1 209 if((ilarger-ismaller)<2 .or. (llarger-lsmaller)<2*tol12)exit 210 itrial=(ilarger+ismaller)/2 ; ltrial=lkpg1_sorted(itrial) 211 if((ltrial-lk2)>tol12)then 212 ilarger=itrial ; llarger=ltrial 213 else if((ltrial-lk2)<-tol12)then 214 ismaller=itrial ; lsmaller=ltrial 215 else 216 ismaller=itrial ; lsmaller=ltrial 217 ilarger=itrial ; llarger=ltrial 218 end if 219 end do 220 itrial=ismaller 221 if(abs(llarger-lk2)<abs(lsmaller-lk2)-tol12)itrial=ilarger 222 if(itrial==0)itrial=ilarger 223 ismaller=itrial ; ilarger=itrial 224 ! write(std_out,*)' listkk : starting search at itrial=',itrial 225 226 dksqmn=huge(one) 227 228 ! The ii index is dummy. This avoids an infinite loop. 229 do ii=1,l3*nkpt1 230 ! do ikpt1=1,nkpt1 231 232 ! If the difference in length between the trial vector and the target vector is bigger 233 ! than the already achieved distance, the search is finished ... 234 ldiff=abs(lkpg1_sorted(itrial)-lk2) 235 236 237 ! DEBUG 238 ! write(std_out,*)' listkk : ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,dksqmn=',ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,dksqmn 239 ! ENDDEBUG 240 if(ldiff**2>dksqmn+tol8)exit 241 242 ! If this k-point has already been examined in a previous batch, skip it 243 ! First, compute the minimum of the difference of length of the sets of associated vectors thanks to Umklapp vectors 244 ! with the target vector 245 ikpt1=(isort(itrial)-1)/l3+1 246 min_l=minval(abs(lkpg1((ikpt1-1)*l3+1:(ikpt1-1)*l3+l3)-lk2)) 247 ! Then compare with the current ldiff 248 249 ! DEBUG 250 ! write(std_out,*)' listkk : ikpt1,min_l,ldiff=',ikpt1,min_l,ldiff 251 ! ENDDEBUG 252 253 if(min_l > ldiff-tol12)then 254 255 ! Now, will examine the trial vector, and the symmetric ones 256 !MG FIXME: 257 ! Here there's a possible problem with the order of symmetries because 258 ! in symkpt, time-reversal is the innermost loop. This can create inconsistencies in the symmetry tables. 259 do itimrev=0,timrev_used 260 do isym=1,nsym_used 261 262 ! Select magnetic characteristic of symmetries 263 if(isppol==1 .and. symafm(isym)==-1)cycle 264 if(isppol==2 .and. symafm(isym)==1)cycle 265 266 ! Compute symmetric point to kpt1 267 if(usesym==1)then 268 ! original code only used transpose(symrel) 269 ! kpt1a(:)=symrel(1,:,isym)*kptns1(1,ikpt1)+& 270 ! & symrel(2,:,isym)*kptns1(2,ikpt1)+& 271 ! & symrel(3,:,isym)*kptns1(3,ikpt1) 272 if (present(use_symrec)) then 273 if (use_symrec) then 274 kpt1a(:) = MATMUL(symmat(:,:,isym),kptns1(:,ikpt1)) 275 else 276 kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1)) 277 end if 278 else 279 kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1)) 280 end if 281 kpt1a(:)=(1-2*itimrev)*kpt1a(:) 282 else 283 kpt1a(:)=kptns1(:,ikpt1) 284 end if 285 286 ! Compute difference with respect to kpt2, modulo a lattice vector 287 dk(:)=kptns2(:,ikpt2)-kpt1a(:) 288 if(usesym==1)then 289 ! The tolerance insure similar behaviour on different platforms 290 ! XG120418 : Actually, *assumes* that the closest point will have reduced 291 ! coordinates differing by less than 1/2 . There might be elongated 292 ! cells where this is not correct ... 293 dkint(:)=nint(dk(:)+tol12) 294 dk(:)=dk(:)-dkint(:) 295 else 296 dkint(:)=0 297 end if 298 299 ! Compute norm of the difference vector, and update kpt1 if better. 300 dksq=gmet(1,1)*dk(1)**2+gmet(2,2)*dk(2)**2+& 301 & gmet(3,3)*dk(3)**2+two*(gmet(2,1)*dk(2)*dk(1)+& 302 & gmet(3,2)*dk(3)*dk(2)+gmet(3,1)*dk(3)*dk(1)) 303 304 if (dksq<dksqmn+tol8) then 305 306 ! If exactly the right point (without using symmetries neither umklapp vector), will exit the search 307 ! Note that in this condition, each coordinate is tested separately, without squaring. So, it is a much stronger 308 ! condition than dksqmn<tol12 309 if(sum(abs(kptns2(:,ikpt2)-kptns1(:,ikpt1)))<3*tol12)then 310 ikpt2_done=1 311 end if 312 313 ! Update in three cases : either if succeeded to have exactly the vector, or the distance is better, 314 ! or the distance is only slightly worsened so select the lowest itimrev, isym or ikpt1, in order to respect previous ordering 315 if( ikpt2_done==1 .or. & 316 & dksq+tol12<dksqmn .or. & 317 & ( abs(dksq-dksqmn)<tol12 .and. & 318 & ((itimrev<jtime) .or. & 319 & (itimrev==jtime .and. isym<jsym) .or. & 320 & (itimrev==jtime .and. isym==jsym .and. ikpt1<jkpt1))))then 321 322 dksqmn=dksq 323 jkpt1=ikpt1 324 jsym=isym 325 jtime=itimrev 326 jdkint(:)=dkint(:) 327 328 ! DEBUG 329 ! write(std_out,*)' ikpt1,ikpt2=',ikpt1,ikpt2 330 ! write(std_out,*)' timrev_used=',timrev_used 331 ! write(std_out,*)' Succeeded to lower dskmn,ikpt2_done=',dksqmn,ikpt2_done 332 ! write(std_out,*)' ikpt1,isym,dkint(:),itimrev=',ikpt1,isym,dkint(:),itimrev 333 ! ka(:)=kpt1a(:)+dkint(:) 334 ! write(std_out,*)' k1=',kpt1a(:) 335 ! write(std_out,*)' dkint=',dkint(:) 336 ! write(std_out,*)' Actual k1=',ka(:) 337 ! write(std_out,*)' k2=',kptns2(:,ikpt2) 338 ! kasq=gmet(1,1)*ka(1)**2+gmet(2,2)*ka(2)**2+& 339 ! & gmet(3,3)*ka(3)**2+two*(gmet(2,1)*ka(2)*ka(1)+& 340 ! & gmet(3,2)*ka(3)*ka(2)+gmet(3,1)*ka(3)*ka(1)) 341 ! write(std_out,*)' Actual k1sq=',kasq 342 ! ENDDEBUG 343 end if 344 345 end if 346 if(ikpt2_done==1)exit 347 348 end do ! isym 349 if(ikpt2_done==1)exit 350 351 end do ! itimrev 352 if(ikpt2_done==1)exit 353 354 end if 355 356 ! Update the interval that has been explored 357 if(itrial<ismaller)ismaller=itrial 358 if(itrial>ilarger)ilarger=itrial 359 360 ! Select the next index to be tried (preferably the smaller indices, but this is a bit arbitrary). 361 362 ! DEBUG 363 ! write(std_out,*)' before choosing the next index :' 364 ! write(std_out,*)' ismaller,itrial,ilarger=',ismaller,itrial,ilarger 365 ! write(std_out,*)' lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1)=',lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1) 366 ! ENDDEBUG 367 if(ismaller>1 .and. ilarger<l3*nkpt1)then 368 if(abs(lkpg1_sorted(ismaller-1)-lk2)<abs(lkpg1_sorted(ilarger+1)-lk2)+tol12)then 369 itrial=ismaller-1 370 else 371 itrial=ilarger+1 372 end if 373 end if 374 if(ismaller==1 .and. ilarger<l3*nkpt1)itrial=ilarger+1 375 if(ismaller>1 .and. ilarger==l3*nkpt1)itrial=ismaller-1 376 ! if(ismaller==1 .and. ilarger==l3*nkpt1), we are done with the loop ! 377 378 end do ! ikpt1 379 380 indkk(ikpt2+(isppol-1)*nkpt2,1)=jkpt1 381 indkk(ikpt2+(isppol-1)*nkpt2,2)=jsym 382 indkk(ikpt2+(isppol-1)*nkpt2,3:5)=jdkint(:) 383 indkk(ikpt2+(isppol-1)*nkpt2,6)=jtime 384 dksqmax=max(dksqmax,dksqmn) 385 386 if(dksqmn<-tol12)then 387 write(message, '(a,es16.6)' )' The minimum square of dk has negative norm: dksqmn=',dksqmn 388 MSG_BUG(message) 389 end if 390 391 ! DEBUG 392 ! write(std_out,'(a,i6,i2,2x,i6,5i3,es24.14)' )' listkk: ikpt2,isppol,indkk(ikpt2+(isppol-1)*nkpt2,:)=',ikpt2,isppol,indkk(ikpt2+(isppol-1)*nkpt2,:),dksqmn 393 ! if(nkpt1==17)stop 394 ! ENDDEBUG 395 396 end do ! ikpt2 397 end do ! isppol 398 399 ABI_DEALLOCATE(isort) 400 ABI_DEALLOCATE(lkpg1) 401 ABI_DEALLOCATE(lkpg1_sorted) 402 403 call timab(1021,2,tsec) 404 405 end subroutine listkk