Branch data 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 : 12332 : 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 : 12332 : integer,allocatable :: list(:),bz2ibz_idx(:)
115 : : real(dp) :: gmetkpt(3),ksym(3)
116 : 12332 : 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 [ - + ]: 12332 : 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 : 12332 : identi = 1
137 : 12332 : tident = -1
138 [ + + ]: 12332 : if (nsym/=1) then
139 [ + - ]: 9310 : do isym=1,nsym
140 : : tident=1
141 [ + + ]: 37240 : do jj=1,3
142 [ - + ]: 27930 : if(symrec(jj,jj,isym)/=1)tident=0
143 [ + + ]: 121030 : do ii=1,3
144 [ + + + - ]: 195510 : if( ii/=jj .and. symrec(ii,jj,isym)/=0)tident=0
145 : : end do
146 : : end do
147 [ - + ]: 9310 : 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 [ - + ]: 9310 : 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 [ + + ]: 5201546 : do ikpt=1,nkbz
158 : 5201546 : wtk_folded(ikpt)=wtk(ikpt)
159 : : end do
160 : :
161 : : ! Initialize bz2ibz_smap
162 [ + + + + ]: 36336830 : bz2ibz_smap = 0
163 [ + + ]: 5201546 : do ikpt=1,nkbz
164 : 5189214 : bz2ibz_smap(1, ikpt) = ikpt
165 : 5201546 : 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 [ + + + + : 12332 : 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 [ - + + - : 29430 : ABI_MALLOC(length2,(nkbz))
+ - + - -
+ - + ]
177 : :
178 [ + + ]: 5070685 : do ikpt=1,nkbz
179 [ + + ]: 163652799 : do isym=1,nsym
180 [ - + + + ]: 478619537 : do itim=1,(1-2*timrev),-2
181 : : ! Get the symmetric of the vector
182 [ + + ]: 1259906192 : do ii=1,3
183 : : ksym(ii)=itim*( kbz(1,ikpt)*symrec(ii,1,isym)&
184 : : +kbz(2,ikpt)*symrec(ii,2,isym)&
185 : 944929644 : +kbz(3,ikpt)*symrec(ii,3,isym) )
186 : 1259906192 : ksym(ii)=ksym(ii)-anint(ksym(ii)+tol8*half)
187 : : end do
188 [ + + ]: 1259906192 : gmetkpt(:)=gmet(:,1)*ksym(1)+gmet(:,2)*ksym(2)+gmet(:,3)*ksym(3)
189 : 314976548 : length2trial=ksym(1)*gmetkpt(1)+ksym(2)*gmetkpt(2)+ksym(3)*gmetkpt(3)
190 [ + + ]: 473558662 : if(isym==1 .and. itim==1)then
191 : 5060875 : length2(ikpt)=length2trial
192 : : else
193 [ + + ]: 309915673 : 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 [ - + + - : 29430 : ABI_MALLOC(list,(nkbz))
+ - + - -
+ - + ]
203 [ + - + + : 10131560 : list(:)=(/ (ikpt,ikpt=1,nkbz) /)
+ + ]
204 : 9810 : 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 [ + + + + ]: 9810 : if (chksymbreak == 1 .and. nkbz < 40**3) then
212 : : ikpt_current_length=1
213 : : ! Loop on all k points
214 [ + + ]: 1171007 : do ikpt=1,nkbz
215 : 1166410 : ind_ikpt=list(ikpt)
216 : : ! Keep track of the current length, to avoid doing needless comparisons
217 [ + + ]: 1166410 : if(length2(ikpt)-length2(ikpt_current_length)>tol8) ikpt_current_length=ikpt
218 : :
219 [ + + ]: 41181719 : do isym=1,nsym
220 [ - + + + ]: 120807026 : do itim=1,(1-2*timrev),-2
221 [ + + ]: 119640616 : if(isym/=identi .or. itim/=1 )then
222 : : ! Get the symmetric of the vector
223 [ + + ]: 313853976 : 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 : 313853976 : +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 : 2656599368 : quit=0
231 [ + - ]: 2656599368 : 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 [ - + ]: 2656599368 : if(length2(ikpt2)-length2(ikpt)>tol8)exit
235 : : ! Ordered index
236 : 2656599368 : ind_ikpt2=list(ikpt2)
237 : 2656599368 : difk1= ksym(1)-kbz(1,ind_ikpt2)
238 : 2656599368 : reduce1=difk1-anint(difk1)
239 : 2656599368 : difk2= ksym(2)-kbz(2,ind_ikpt2)
240 : 2656599368 : reduce2=difk2-anint(difk2)
241 : 2656599368 : difk3= ksym(3)-kbz(3,ind_ikpt2)
242 : 2656599368 : reduce3=difk3-anint(difk3)
243 [ + + ]: 2656599368 : 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 [ - + ]: 78463494 : 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 [ + + ]: 5060875 : do ikpt=1,nkbz-1
270 : : ! Ordered index
271 : 5051065 : 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 [ + + ]: 5051065 : if (wtk_folded(ind_ikpt) < tol16) cycle
276 : :
277 : : ! Loop on the remaining k-points
278 [ + + ]: 14293215 : 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 [ + + ]: 14267565 : if (length2(ikpt2) - length2(ikpt) > tol8) exit
284 : :
285 : : ! Ordered index
286 : 13810038 : ind_ikpt2=list(ikpt2)
287 : :
288 : : ! If the second vector is already empty, no interest to treat it
289 [ + + ]: 13810038 : if (wtk_folded(ind_ikpt2) < tol16) cycle
290 : :
291 : 10122669 : 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 [ + + ]: 236559693 : do isym=1,nsym
296 [ - + + + ]: 668939267 : do itim=1,(1-2*timrev),-2
297 [ + + + + ]: 668939267 : if (isym/=identi .or. itim/=1) then
298 : : ! Get the symmetric of the vector
299 [ + + ]: 1749722556 : 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 : 1749722556 : +kbz(3,ind_ikpt)*symrec(ii,3,isym) )
303 : : end do
304 : :
305 : : ! The do-loop was expanded to speed up the execution
306 : 437430639 : difk= ksym(1)-kbz(1,ind_ikpt2)
307 : 437430639 : reduce=difk-anint(difk)
308 [ + + ]: 437430639 : if (abs(reduce)>tol8) cycle
309 : 29410380 : difk= ksym(2)-kbz(2,ind_ikpt2)
310 : 29410380 : reduce=difk-anint(difk)
311 [ + + ]: 29410380 : if (abs(reduce)>tol8) cycle
312 : 7141508 : difk= ksym(3)-kbz(3,ind_ikpt2)
313 : 7141508 : reduce=difk-anint(difk)
314 [ + + ]: 7141508 : 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 : 4585793 : wtk_folded(ind_ikpt) = wtk_folded(ind_ikpt) + wtk_folded(ind_ikpt2)
319 : 4585793 : 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 : 4585793 : bz2ibz_smap(1, ind_ikpt2) = ind_ikpt
325 : 4585793 : bz2ibz_smap(2, ind_ikpt2) = isym
326 : : ! Compute difference with respect to kpt2, modulo a lattice vector
327 : : ! Sk1 + G0 = k2
328 [ + + ]: 18343172 : bz2ibz_smap(3:5, ind_ikpt2) = nint(-ksym(:) + kbz(:, ind_ikpt2) + tol12)
329 [ + + ]: 4585793 : ii = 0; if (itim == -1) ii = 1
330 : 4585793 : 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 : 9224245 : if (quit == 1) exit
337 : : end do ! itim
338 : :
339 : : end do ! ikpt2
340 : : end do ! ikpt
341 : :
342 [ - + ]: 9810 : ABI_FREE(length2)
343 [ - + ]: 9810 : 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 [ - + + - : 36996 : ABI_MALLOC(bz2ibz_idx, (nkbz))
+ - + - -
+ - + ]
349 [ + + ]: 5201546 : bz2ibz_idx = 0
350 : 12332 : nkibz = 0
351 [ + + ]: 5201546 : do ikpt=1,nkbz
352 [ + + ]: 5201546 : if (wtk_folded(ikpt) > tol8) then
353 : 603421 : nkibz = nkibz+1
354 : 603421 : ibz2bz(nkibz) = ikpt
355 : 603421 : 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 [ + + ]: 5201546 : do ikpt=1,nkbz
362 : 5189214 : ind_ikpt = bz2ibz_idx(bz2ibz_smap(1, ikpt))
363 [ + - ]: 5201546 : if (ind_ikpt /= 0) then
364 : 5189214 : bz2ibz_smap(1, ikpt) = ind_ikpt
365 : : else
366 : 0 : ierr = ierr + 1
367 : : end if
368 : : end do
369 [ - + ]: 12332 : ABI_CHECK(ierr == 0, "Error while remapping bz2ibz_smap array")
370 : :
371 [ - + ]: 12332 : ABI_FREE(bz2ibz_idx)
372 : :
373 [ + + ]: 12332 : if(iout/=0)then
374 [ + + ]: 2614 : if(nkbz/=nkibz)then
375 : : write(msg, '(a,a,a,i6,a)' )&
376 : 636 : ' symkpt : the number of k-points, thanks to the symmetries,',ch10,' is reduced to',nkibz,' .'
377 : 636 : call wrtout(iout,msg)
378 [ + + ]: 636 : if(iout/=std_out) call wrtout(std_out,msg)
379 : :
380 : 636 : 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 : 1978 : write(msg, '(a)' )' symkpt : not enough symmetry to change the number of k points.'
399 : 1978 : call wrtout(iout,msg)
400 [ + + ]: 1978 : if (iout/=std_out) call wrtout(std_out,msg)
401 : : end if
402 : : end if
403 : :
404 : 12332 : 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 : : !!***
|
