Line data Source code
1 : !!****m* ABINIT/m_symkpt
2 : !! NAME
3 : !! m_symkpt
4 : !!
5 : !! COPYRIGHT
6 : !! Copyright (C) 1999-2026 ABINIT group (XG,LSI,HM)
7 : !! This file is distributed under the terms of the
8 : !! GNU General Public License, see ~abinit/COPYING
9 : !! or http://www.gnu.org/copyleft/gpl.txt .
10 : !!
11 : !! TODO
12 : !! Move it to m_kpts
13 : !!
14 : !! SOURCE
15 :
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
37 : !!***
38 :
39 : public :: symkpt
40 : public :: symkpt_new
41 : public :: mapkptsets
42 : !!***
43 :
44 : contains
45 : !!***
46 :
47 : !!****f* ABINIT/symkpt
48 : !! NAME
49 : !! symkpt
50 : !!
51 : !! FUNCTION
52 : !! Determines the weights of the k-points for sampling the Brillouin Zone, starting from a first set
53 : !! of weights wtk, and folding it to a new set, by taking into account the symmetries described
54 : !! by symrec, and eventually the time-reversal symmetry.
55 : !! Also compute the number of k points in the reduced set
56 : !! This routine is also used for sampling the q vectors in the Brillouin zone for the computation
57 : !! of thermodynamical properties (from the routine thm9).
58 : !!
59 : !! INPUTS
60 : !! chksymbreak= if 1, will check whether the k point grid is symmetric, and stop if not.
61 : !! gmet(3,3)=reciprocal space metric (bohr**-2).
62 : !! iout=if non-zero, output the new number of kpoints on unit iout
63 : !! kbz(3,nkbz)= k vectors in the BZ.
64 : !! nkbz = number of k-points whose weights are wtk
65 : !! nsym=number of space group symmetries
66 : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
67 : !! timrev: if 1, the time reversal operation has to be taken into account
68 : !! if 0, no time reversal symmetry.
69 : !! wtk(nkbz)=weight assigned to each k point.
70 : !! comm=MPI communicator.
71 : !!
72 : !! OUTPUT
73 : !! ibz2bz(nkbz)=non-symmetrized indices of the k-points a.k.a. ibz2bz mapping
74 : !! The correspondence beween the iq_ibz point in IBZ and the iq_bz point in the full BZ is obtained via:
75 : !!
76 : !! do ik_ibz=1,nkibz
77 : !! ik_bz = ibz2bz(ik_ibz)
78 : !! end do
79 : !!
80 : !! nkibz = number of k-points in the irreducible set
81 : !! wtk_folded(nkbz)=weight assigned to each k point, taking into account the symmetries
82 : !! bz2ibz_smap(nkbz, 6)= Mapping BZ --> IBZ.
83 : !!
84 : !! NOTES
85 : !! The decomposition of the symmetry group in its primitives might speed up the execution.
86 : !! The output variables are stored only in the range 1:nkbz
87 : !!
88 : !! TODO
89 : !! Bad scaling wrt nkbz. Should try to MPI parallelize or implement more efficient algorithm
90 : !!
91 : !! SOURCE
92 :
93 12248 : 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 12248 : integer,allocatable :: list(:),bz2ibz_idx(:)
115 : real(dp) :: gmetkpt(3),ksym(3)
116 12248 : 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 12248 : if (timrev/=1 .and. timrev/=0) then
132 0 : 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 12248 : identi = 1
137 12248 : tident = -1
138 12248 : if (nsym/=1) then
139 9258 : do isym=1,nsym
140 : tident=1
141 37032 : do jj=1,3
142 27774 : if(symrec(jj,jj,isym)/=1)tident=0
143 120354 : do ii=1,3
144 194418 : if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
145 : end do
146 : end do
147 9258 : 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 9258 : 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 5779878 : do ikpt=1,nkbz
158 5779878 : wtk_folded(ikpt)=wtk(ikpt)
159 : end do
160 :
161 : ! Initialize bz2ibz_smap
162 40385658 : bz2ibz_smap = 0
163 5779878 : do ikpt=1,nkbz
164 5767630 : bz2ibz_smap(1, ikpt) = ikpt
165 5779878 : 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 12248 : 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 29262 : ABI_MALLOC(length2,(nkbz))
177 :
178 5649157 : do ikpt=1,nkbz
179 188548483 : do isym=1,nsym
180 552153397 : do itim=1,(1-2*timrev),-2
181 : ! Get the symmetric of the vector
182 1454458672 : do ii=1,3
183 : ksym(ii)=itim*( kbz(1,ikpt)*symrec(ii,1,isym)&
184 : +kbz(2,ikpt)*symrec(ii,2,isym)&
185 1090844004 : +kbz(3,ikpt)*symrec(ii,3,isym) )
186 1454458672 : ksym(ii)=ksym(ii)-anint(ksym(ii)+tol8*half)
187 : end do
188 1454458672 : gmetkpt(:)=gmet(:,1)*ksym(1)+gmet(:,2)*ksym(2)+gmet(:,3)*ksym(3)
189 363614668 : length2trial=ksym(1)*gmetkpt(1)+ksym(2)*gmetkpt(2)+ksym(3)*gmetkpt(3)
190 546513994 : if(isym==1 .and. itim==1)then
191 5639403 : length2(ikpt)=length2trial
192 : else
193 357975265 : 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 29262 : ABI_MALLOC(list,(nkbz))
203 11288560 : list(:)=(/ (ikpt,ikpt=1,nkbz) /)
204 9754 : 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 9754 : if (chksymbreak == 1 .and. nkbz < 40**3) then
212 : ikpt_current_length=1
213 : ! Loop on all k points
214 1170175 : do ikpt=1,nkbz
215 1165610 : ind_ikpt=list(ikpt)
216 : ! Keep track of the current length, to avoid doing needless comparisons
217 1165610 : if(length2(ikpt)-length2(ikpt_current_length)>tol8) ikpt_current_length=ikpt
218 :
219 41149959 : do isym=1,nsym
220 120713442 : do itim=1,(1-2*timrev),-2
221 119547832 : if(isym/=identi .or. itim/=1 )then
222 : ! Get the symmetric of the vector
223 313609752 : 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 313609752 : +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 2656163264 : quit=0
231 2656163264 : 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 2656163264 : if(length2(ikpt2)-length2(ikpt)>tol8)exit
235 : ! Ordered index
236 2656163264 : ind_ikpt2=list(ikpt2)
237 2656163264 : difk1= ksym(1)-kbz(1,ind_ikpt2)
238 2656163264 : reduce1=difk1-anint(difk1)
239 2656163264 : difk2= ksym(2)-kbz(2,ind_ikpt2)
240 2656163264 : reduce2=difk2-anint(difk2)
241 2656163264 : difk3= ksym(3)-kbz(3,ind_ikpt2)
242 2656163264 : reduce3=difk3-anint(difk3)
243 2656163264 : 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 78402438 : if (quit == 0) then
250 : write(msg,'(3a,i0,2a,9(i0,1x),2a,i0,1a,3es16.6,6a)' )&
251 0 : 'Chksymbreak = 1. It has been observed that the k point grid is not symmetric:',ch10,&
252 0 : 'for the symmetry number: ',isym,ch10,&
253 0 : 'with symrec= ',symrec(1:3,1:3,isym),ch10,&
254 0 : 'the symmetric of the k point number: ',ind_ikpt2,' with components: ', kbz(1:3,ind_ikpt2),ch10,&
255 0 : 'does not belong to the k point grid.',ch10,&
256 0 : 'Read the description of the input variable chksymbreak,',ch10,&
257 0 : 'You might switch it to zero, or change your k point grid to one that is symmetric.'
258 0 : 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 5639403 : do ikpt=1,nkbz-1
270 : ! Ordered index
271 5629649 : 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 5629649 : if (wtk_folded(ind_ikpt) < tol16) cycle
276 :
277 : ! Loop on the remaining k-points
278 15947304 : 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 15921790 : if (length2(ikpt2) - length2(ikpt) > tol8) exit
284 :
285 : ! Ordered index
286 15392886 : ind_ikpt2=list(ikpt2)
287 :
288 : ! If the second vector is already empty, no interest to treat it
289 15392886 : if (wtk_folded(ind_ikpt2) < tol16) cycle
290 :
291 11224891 : 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 284919890 : do isym=1,nsym
296 809776580 : do itim=1,(1-2*timrev),-2
297 809776580 : if (isym/=identi .or. itim/=1) then
298 : ! Get the symmetric of the vector
299 2121945356 : 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 2121945356 : +kbz(3,ind_ikpt)*symrec(ii,3,isym) )
303 : end do
304 :
305 : ! The do-loop was expanded to speed up the execution
306 530486339 : difk= ksym(1)-kbz(1,ind_ikpt2)
307 530486339 : reduce=difk-anint(difk)
308 530486339 : if (abs(reduce)>tol8) cycle
309 34504276 : difk= ksym(2)-kbz(2,ind_ikpt2)
310 34504276 : reduce=difk-anint(difk)
311 34504276 : if (abs(reduce)>tol8) cycle
312 8229101 : difk= ksym(3)-kbz(3,ind_ikpt2)
313 8229101 : reduce=difk-anint(difk)
314 8229101 : 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 5093028 : wtk_folded(ind_ikpt) = wtk_folded(ind_ikpt) + wtk_folded(ind_ikpt2)
319 5093028 : 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 5093028 : bz2ibz_smap(1, ind_ikpt2) = ind_ikpt
325 5093028 : bz2ibz_smap(2, ind_ikpt2) = isym
326 : ! Compute difference with respect to kpt2, modulo a lattice vector
327 : ! Sk1 + G0 = k2
328 20372112 : bz2ibz_smap(3:5, ind_ikpt2) = nint(-ksym(:) + kbz(:, ind_ikpt2) + tol12)
329 5093028 : ii = 0; if (itim == -1) ii = 1
330 5093028 : 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 10299858 : if (quit == 1) exit
337 : end do ! itim
338 :
339 : end do ! ikpt2
340 : end do ! ikpt
341 :
342 9754 : ABI_FREE(length2)
343 9754 : 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 36744 : ABI_MALLOC(bz2ibz_idx, (nkbz))
349 5779878 : bz2ibz_idx = 0
350 12248 : nkibz = 0
351 5779878 : do ikpt=1,nkbz
352 5779878 : if (wtk_folded(ikpt) > tol8) then
353 674602 : nkibz = nkibz+1
354 674602 : ibz2bz(nkibz) = ikpt
355 674602 : 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 5779878 : do ikpt=1,nkbz
362 5767630 : ind_ikpt = bz2ibz_idx(bz2ibz_smap(1, ikpt))
363 5779878 : if (ind_ikpt /= 0) then
364 5767630 : bz2ibz_smap(1, ikpt) = ind_ikpt
365 : else
366 0 : ierr = ierr + 1
367 : end if
368 : end do
369 12248 : ABI_CHECK(ierr == 0, "Error while remapping bz2ibz_smap array")
370 :
371 12248 : ABI_FREE(bz2ibz_idx)
372 :
373 12248 : if(iout/=0)then
374 2598 : if(nkbz/=nkibz)then
375 : write(msg, '(a,a,a,i6,a)' )&
376 632 : ' symkpt : the number of k-points, thanks to the symmetries,',ch10,' is reduced to',nkibz,' .'
377 632 : call wrtout(iout,msg)
378 632 : if(iout/=std_out) call wrtout(std_out,msg)
379 :
380 632 : 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 1966 : write(msg, '(a)' )' symkpt : not enough symmetry to change the number of k points.'
399 1966 : call wrtout(iout,msg)
400 1966 : if (iout/=std_out) call wrtout(std_out,msg)
401 : end if
402 : end if
403 :
404 12248 : end subroutine symkpt
405 : !!***
406 :
407 : !----------------------------------------------------------------------
408 :
409 : !!****f* ABINIT/symkpt_new
410 : !! NAME
411 : !! symkpt_new
412 : !!
413 : !! FUNCTION
414 : !! Same routine as above but with an algorithm with better scalling than before.
415 : !! From a few tests it produces the same IBZ as before but avoids computing the lengths and sorting.
416 : !! Instead it uses the krank datatype to map k-points onto each other.
417 : !!
418 : !! INPUTS
419 : !!
420 : !! OUTPUT
421 : !!
422 : !! SOURCE
423 :
424 209 : 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 209 : 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 209 : if (timrev/=1 .and. timrev/=0) then
452 0 : write(msg,'(a,i0)')' timrev should be 0 or 1, while it is equal to ',timrev
453 0 : ABI_BUG(msg)
454 : end if
455 :
456 : ! Find the identity symmetry operation
457 209 : identi = 1
458 209 : tident = -1
459 209 : if (nsym/=1) then
460 209 : do isym=1,nsym
461 : tident=1
462 836 : do jj=1,3
463 627 : if(symrec(jj,jj,isym)/=1)tident=0
464 2717 : do ii=1,3
465 4389 : if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
466 : end do
467 : end do
468 209 : 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 209 : ABI_CHECK(tident == 1, 'Did not find the identity operation')
475 : end if
476 :
477 : ! Initialize
478 820945 : ibz2bz = 0
479 5745361 : bz2ibz_smap = 0
480 820945 : do ikpt=1,nkbz
481 820736 : bz2ibz_smap(1, ikpt) = ikpt
482 820736 : bz2ibz_smap(2, ikpt) = 1
483 820945 : bz2ibz_smap(4, ikpt) = 1 ! We will use this as wtk_folded
484 : end do
485 :
486 : ! Start krank
487 209 : call krank%init(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 209 : 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 209 : if (chksymbreak == 1 .and. nkbz < 40**3) then
498 0 : do ikpt=1,nkbz
499 0 : kpt1 = kbz(:,ikpt)
500 :
501 0 : do isym=1,nsym
502 0 : do itim=0,timrev
503 : ! Skip identity symmetry
504 0 : if (isym==identi .and. itim==0) cycle
505 :
506 : ! Get the symmetric of the vector
507 0 : do ii=1,3
508 : ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
509 : kpt1(2)*symrec(ii,2,isym)+&
510 0 : kpt1(3)*symrec(ii,3,isym) )
511 : end do
512 :
513 : !find this point
514 0 : 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 0 : if (ikpt_found < 0) then
520 : write(msg,'(3a,i4,2a,9i3,2a,i6,1a,3es16.6,6a)' )&
521 0 : 'Chksymbreak=1. It has been observed that the k point grid is not symmetric:',ch10,&
522 0 : 'for the symmetry number: ',isym,ch10,&
523 0 : 'with symrec= ',symrec(1:3,1:3,isym),ch10,&
524 0 : 'the symmetric of the k point number: ',ikpt,' with components: ',kpt1(:),ch10,&
525 0 : 'does not belong to the k point grid.',ch10,&
526 0 : 'Read the description of the input variable chksymbreak,',ch10,&
527 0 : 'You might switch it to zero, or change your k point grid to one that is symmetric.'
528 0 : 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 820945 : ikpt_loop: do ikpt=1,nkbz
537 :
538 820736 : if (bz2ibz_smap(4, ikpt)==0) cycle
539 3282944 : 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 8320081 : do isym=1,nsym
545 15819008 : do itim=0,timrev
546 : ! Skip identity symmetry
547 7499136 : if (isym==identi .and. itim==0) cycle
548 :
549 : ! Get the symmetric of the vector
550 26713600 : do ii=1,3
551 : ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
552 : kpt1(2)*symrec(ii,2,isym)+&
553 26713600 : kpt1(3)*symrec(ii,3,isym) )
554 : end do
555 :
556 : !find this point
557 6678400 : ikpt_found = krank%get_index(ksym)
558 : !if k-point not found just cycle
559 6678400 : 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 6678400 : if (ikpt_found >= ikpt) cycle
564 12969840 : bz2ibz_smap(:3, ikpt) = [ikpt_found, isym, itim]
565 3242460 : bz2ibz_smap(4,ikpt) = bz2ibz_smap(4,ikpt) + bz2ibz_smap(4,ikpt_found)
566 14998272 : 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 209 : nkibz = 0
575 820945 : do ikpt=1,nkbz
576 820736 : ikibz = bz2ibz_smap(1,ikpt)
577 820736 : if (ikibz /= ikpt) cycle
578 214821 : nkibz = nkibz + 1
579 820945 : 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 5745361 : bz2ibz_smap = 0
588 :
589 : ! Now I loop again over the points in the IBZ to find the mapping to the BZ
590 215030 : do ikpt=1,nkibz
591 214821 : ikibz = ibz2bz(ikpt)
592 859284 : kpt1 = kbz(:,ikibz)
593 :
594 : ! HM: Here I invert the itim and isym loop to generate the same mapping as listkk
595 429851 : do itim=0,timrev
596 1443858 : do isym=1,nsym
597 : ! Get the symmetric of the vector
598 4056864 : do ii=1,3
599 : ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
600 : kpt1(2)*symrec(ii,2,isym)+&
601 4056864 : kpt1(3)*symrec(ii,3,isym) )
602 : end do
603 :
604 : ! Find this point
605 1014216 : ikpt_found = krank%get_index(ksym)
606 1014216 : if (ikpt_found < 0) cycle
607 1014216 : if (bz2ibz_smap(1, ikpt_found) /= 0) cycle
608 3282944 : bz2ibz_smap(:3, ikpt_found) = [ikpt, isym, itim]
609 3691245 : 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 209 : call krank%free()
621 :
622 : ! Check whether the mapping was sucessfull
623 820945 : if (any(bz2ibz_smap(1, :) == 0)) then
624 0 : 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 209 : if(iout/=0)then
632 0 : if(nkbz/=nkibz)then
633 : write(msg, '(a,a,a,i6,a)' )&
634 0 : ' symkpt_new : the number of k-points, thanks to the symmetries,',ch10,' is reduced to',nkibz,' .'
635 0 : call wrtout(iout,msg)
636 0 : if(iout/=std_out) call wrtout(std_out,msg)
637 :
638 0 : nkpout=nkibz
639 : else
640 0 : write(msg, '(a)' )' symkpt_new : not enough symmetry to change the number of k points.'
641 0 : call wrtout(iout,msg)
642 0 : if (iout/=std_out) call wrtout(std_out,msg)
643 : end if
644 : end if
645 :
646 209 : end subroutine symkpt_new
647 : !!***
648 :
649 : !----------------------------------------------------------------------
650 :
651 : !!****f* ABINIT/mapkptsets
652 : !! NAME
653 : !! mapkptsets
654 : !!
655 : !! FUNCTION
656 : !! 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
657 : !! typical usage is to find k in a list from disk, to initialize the wfk in memory
658 : !! kin can be overcomplete etc... we just need to find _a_ solution
659 : !! bz2kin_smap follows the listkk convention with ik, isym, g0(1:3), itimrev
660 : !!
661 : !! INPUTS
662 : !!
663 : !! OUTPUT
664 : !! nkirred = number of irreducible k needed from list 1 (in)
665 : !! bz2ibz_smap = mapping of indices in list 2 with their irreducible origin in list 1,
666 : !! symop and timrev needed to transform them
667 : !!
668 : !! SOURCE
669 :
670 0 : subroutine mapkptsets(chksymbreak,gmet,k_in,nk_in,&
671 0 : & 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 0 : 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 0 : if (timrev/=1 .and. timrev/=0) then
700 0 : write(message,'(a,i0)')' timrev should be 0 or 1, while it is equal to ',timrev
701 0 : ABI_BUG(message)
702 : end if
703 :
704 : ! Find the identity symmetry operation
705 0 : identi = 1
706 0 : tident = -1
707 0 : if (nsym/=1) then
708 0 : do isym=1,nsym
709 : tident=1
710 0 : do jj=1,3
711 0 : if(symrec(jj,jj,isym)/=1)tident=0
712 0 : do ii=1,3
713 0 : if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
714 : end do
715 : end do
716 0 : if(tident==1)then
717 : identi=isym
718 : exit
719 : end if
720 : end do
721 0 : ABI_CHECK(tident == 1, 'Did not find the identity operation')
722 : end if
723 :
724 : ! Initialize
725 0 : bz2kin_smap = 0
726 0 : do ikpt=1,nkbz
727 0 : bz2kin_smap(ikpt, 1) = ikpt
728 0 : bz2kin_smap(ikpt, 2) = 1
729 0 : bz2kin_smap(ikpt, 3) = 1 ! We will use this as wtk_folded
730 : end do
731 :
732 : ! Start krank
733 0 : call krank%init(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 0 : 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 0 : if (chksymbreak == 1 .and. nkbz < 40**3) then
744 0 : do ikpt=1,nkbz
745 0 : kpt1 = kbz(:,ikpt)
746 :
747 0 : do isym=1,nsym
748 0 : do itim=0,timrev
749 : ! Skip identity symmetry
750 0 : if (isym==identi .and. itim==0) cycle
751 :
752 : ! Get the symmetric of the vector
753 0 : do ii=1,3
754 : ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
755 : kpt1(2)*symrec(ii,2,isym)+&
756 0 : kpt1(3)*symrec(ii,3,isym) )
757 : end do
758 :
759 : !find this point
760 0 : 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 0 : if (ikpt_found < 0) then
766 : write(message,'(3a,i4,2a,9i3,2a,i6,1a,3es16.6,6a)' )&
767 0 : 'Chksymbreak=1. It has been observed that the k point grid is not symmetric:',ch10,&
768 0 : 'for the symmetry number: ',isym,ch10,&
769 0 : 'with symrec= ',symrec(1:3,1:3,isym),ch10,&
770 0 : 'the symmetric of the k point number: ',ikpt,' with components: ',kpt1(:),ch10,&
771 0 : 'does not belong to the k point grid.',ch10,&
772 0 : 'Read the description of the input variable chksymbreak,',ch10,&
773 0 : 'You might switch it to zero, or change your k point grid to one that is symmetric.'
774 0 : 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 0 : bz2kin_smap = 0
785 :
786 : ! HM: Here I invert the itim and isym loop to generate the same mapping as listkk
787 0 : do itim=0,timrev
788 0 : 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 0 : do ik_in=1,nk_in
794 0 : kpt1 = k_in(:,ik_in)
795 :
796 : ! Get the symmetric of the vector
797 0 : do ii=1,3
798 : ksym(ii)=(1-2*itim)*( kpt1(1)*symrec(ii,1,isym)+&
799 : kpt1(2)*symrec(ii,2,isym)+&
800 0 : kpt1(3)*symrec(ii,3,isym) )
801 : end do
802 :
803 : !find this point in the main set 2
804 0 : ikpt_found = krank%get_index(ksym)
805 :
806 0 : if (ikpt_found < 0) cycle
807 : ! if we already filled it, ignore new symmetric pre-image
808 0 : if (bz2kin_smap(ikpt_found, 1) /= 0) cycle
809 :
810 0 : bz2kin_smap(ikpt_found, 1) = ik_in
811 0 : bz2kin_smap(ikpt_found, 2) = isym
812 0 : bz2kin_smap(ikpt_found, 3:5) = nint(kbz(:,ikpt_found)-ksym)
813 0 : 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 0 : nkirred = 0
821 0 : do ik_in=1,nk_in
822 : ! did I end up using this ik_in?
823 0 : if (any(bz2kin_smap(:,1) == ik_in)) then
824 0 : nkirred = nkirred + 1
825 : end if
826 : end do
827 :
828 : ! check for redundant k in the kbz set
829 0 : do ikpt=1, nkbz
830 0 : if (bz2kin_smap(ikpt,1) /= 0) cycle
831 : ! find index according to krank
832 0 : ikpt_found = krank%get_index(kbz(:,ikpt))
833 : ! if I am not my own hash image, associate the smap data from my image
834 0 : if (ikpt_found /= ikpt) then
835 0 : bz2kin_smap(ikpt,:) = bz2kin_smap(ikpt_found,:)
836 : end if
837 : end do
838 :
839 0 : 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 0 : 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 0 : 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 0 : end subroutine mapkptsets
860 : !!***
861 :
862 : end module m_symkpt
863 : !!***
|
