Branch data Line data Source code
1 : : !!****m* ABINIT/m_krank
2 : : !! NAME
3 : : !! m_krank
4 : : !!
5 : : !! FUNCTION
6 : : !! This module deals with rank objects for hashing k-point vector lists
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2010-2026 ABINIT group (MVer, HM, MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public Licence, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_krank
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_errors
28 : :
29 : : use m_fstrings, only : itoa, sjoin
30 : :
31 : : implicit none
32 : :
33 : : private
34 : : !!***
35 : :
36 : : !!****t* m_krank/krank_t
37 : : !! NAME
38 : : !! krank_t
39 : : !!
40 : : !! FUNCTION
41 : : !! structure to contain a rank/inverse rank pair of arrays, with dimensions
42 : : !!
43 : : !! SOURCE
44 : :
45 : : type,public :: krank_t
46 : :
47 : : integer :: max_linear_density = -1
48 : :
49 : : integer :: min_rank = -1
50 : :
51 : : integer :: max_rank = -1
52 : :
53 : : integer :: npoints = -1
54 : :
55 : : logical :: time_reversal
56 : :
57 : : logical :: kpts_owns_memory = .False.
58 : :
59 : : integer,allocatable :: invrank(:)
60 : :
61 : : real(dp),contiguous, pointer :: kpts(:,:)
62 : : ! Reference to input k-points or copy of the array depending on kpts_owns_memory
63 : :
64 : : ! Internal tables used by krank_get_mapping
65 : : integer,allocatable :: rank2ikpt_(:), rank2symtime_(:)
66 : :
67 : : contains
68 : :
69 : : procedure :: init => krank_init
70 : : ! Sets up the kpt ranks for comparing kpts
71 : :
72 : : procedure :: from_kptrlatt => krank_from_kptrlatt
73 : : ! Initialize object from kptrlatt
74 : :
75 : : procedure :: get_rank
76 : : ! Calculates the rank for one kpt
77 : :
78 : : procedure :: get_index => krank_get_index
79 : : ! Return the index of the k-point `kpt` in the initial set. -1 if not found.
80 : :
81 : : procedure :: copy => krank_copy
82 : : ! Deep copy of the object
83 : :
84 : : procedure :: free => krank_free
85 : : ! Free memory
86 : :
87 : : procedure :: print => krank_print
88 : : ! Prints the arrays and dimensions of a krank_t structure
89 : :
90 : : procedure :: get_mapping => krank_get_mapping
91 : : ! Use symmetries to map input kptn2 to the list of k-points used to generate krank_t.
92 : : ! Similar to listkk but, unlike listkk, this algo does not try to minimize the distance.
93 : : ! Must faster than listkk for dense meshes
94 : :
95 : : end type krank_t
96 : :
97 : : public :: get_ibz2bz ! Return array with the index of the IBZ wave vectors in the BZ.
98 : : public :: star_from_ibz_idx ! Return array with the indices of the star of the ik_ibz wavevector in the IBZ.
99 : : !!***
100 : :
101 : : contains
102 : :
103 : :
104 : : !!****f* m_krank/krank_from_kptrlatt
105 : : !! NAME
106 : : !! krank_from_kptrlatt
107 : : !!
108 : : !! FUNCTION
109 : : !! This routine sets up the kpt ranks for comparing kpts
110 : : !!
111 : : !! INPUTS
112 : : !! npt = number of kpoints (eventually irreducible)
113 : : !! kpt = coordinates of kpoints
114 : : !!
115 : : !! NOTES
116 : : !! By default, the object holds a reference to kpts so do not change/deallocate this array
117 : : !! while using krank.
118 : : !!
119 : : !! SOURCE
120 : :
121 [ + - ]: 478 : subroutine krank_from_kptrlatt(new, nkpt, kpts, kptrlatt, compute_invrank)
122 : :
123 : : !Arguments ------------------------------------
124 : : !scalars
125 : : class(krank_t), intent(out) :: new
126 : : integer,intent(in) :: nkpt
127 : : logical,optional,intent(in) :: compute_invrank
128 : : !arrays
129 : : integer,intent(in) :: kptrlatt(3,3)
130 : : real(dp),target,intent(in) :: kpts(3,nkpt)
131 : :
132 : : !Local variables -------------------------
133 : : !scalars
134 : : integer :: ii, jj, ikpt, max_linear_density, opt=0
135 : : logical :: compute_invrank_
136 : : real(dp) :: min_kpt
137 : : ! *********************************************************************
138 : :
139 : 478 : opt=0
140 [ + + ]: 1912 : do jj=1,3
141 [ + + ]: 6214 : do ii=1,3
142 [ + + - + ]: 4302 : if (ii == jj .and. kptrlatt(ii, ii) == 0) then
143 : 0 : ABI_ERROR("kptrlatt with zero matrix element on the diagonal!")
144 : : end if
145 [ + + - + ]: 5736 : if (ii /= jj .and. kptrlatt(ii, jj) /= 0) then
146 : : !ABI_WARNING("kptrlatt with non-zero off-diagonal matrix elements is not supported")
147 : 0 : opt=1
148 : : end if
149 : : end do
150 : : end do
151 : :
152 [ + - ]: 478 : compute_invrank_ = .True.; if (present(compute_invrank)) compute_invrank_ = compute_invrank
153 : :
154 : 478 : min_kpt = 1
155 [ - + ]: 478 : if (opt == 1) then
156 [ # # ]: 0 : do ikpt=1,nkpt
157 [ # # ]: 0 : do ii=1,3
158 [ # # # # ]: 0 : if (abs(kpts(ii,ikpt)) < min_kpt .and. abs(kpts(ii,ikpt)) /= 0) then
159 : 0 : min_kpt = abs(kpts(ii,ikpt)) ! used as tmp variable
160 : : end if
161 : : end do
162 : : end do
163 [ # # ]: 0 : max_linear_density = ceiling(2/min_kpt)
164 : : else
165 [ + + ]: 1912 : max_linear_density = maxval([kptrlatt(1,1), kptrlatt(2,2), kptrlatt(3,3)])
166 : : end if
167 : :
168 : 478 : call new%init(nkpt, kpts, max_linear_density=max_linear_density, compute_invrank=compute_invrank_)
169 : :
170 : 478 : end subroutine krank_from_kptrlatt
171 : : !!***
172 : :
173 : : !!****f* m_krank/krank_init
174 : : !! NAME
175 : : !! krank_init
176 : : !!
177 : : !! FUNCTION
178 : : !! This routine sets up the kpt ranks for comparing kpts
179 : : !!
180 : : !! INPUTS
181 : : !! npt = number of kpoints (eventually irreducible)
182 : : !! kpt = coordinates of kpoints
183 : : !! time_reversal = true or false to use time reversal symmetry.
184 : : !! Default is true, but only important if nsym and symrec are present
185 : : !!
186 : : !! SOURCE
187 : :
188 [ + - ]: 1113 : subroutine krank_init(new, nkpt, kpts, nsym, symrec, time_reversal, max_linear_density, compute_invrank)
189 : :
190 : : !Arguments ------------------------------------
191 : : !scalars
192 : : class(krank_t),intent(out) :: new
193 : : integer,intent(in) :: nkpt
194 : : integer,intent(in), optional :: nsym
195 : : logical,intent(in), optional :: time_reversal
196 : : integer,optional,intent(in) :: max_linear_density
197 : : logical,optional,intent(in) :: compute_invrank
198 : : !arrays
199 : : real(dp),target,intent(in) :: kpts(3,nkpt)
200 : : integer,intent(in), optional :: symrec(3,3, *)
201 : :
202 : : !Local variables -------------------------
203 : : !scalars
204 : : integer :: ikpt, isym, symkptrank, irank, timrev, itim
205 : : logical :: compute_invrank_
206 : : real(dp) :: smallestlen
207 : : character(len=500) :: msg
208 : : !arrays
209 : : real(dp) :: symkpt(3)
210 : : ! *********************************************************************
211 : :
212 [ + + ]: 1113 : compute_invrank_ = .True.; if (present(compute_invrank)) compute_invrank_ = compute_invrank
213 : :
214 : 1113 : new%kpts => kpts
215 : 1113 : new%kpts_owns_memory = .False.
216 : :
217 [ + + ]: 1113 : if (.not. present(max_linear_density)) then
218 : : ! Find smallest linear length from input kpts
219 : 635 : smallestlen = one
220 [ + + ]: 2358271 : do ikpt=1, nkpt
221 [ + + ]: 2357636 : if (abs(kpts(1,ikpt)) > tol10) smallestlen = min(smallestlen, abs(kpts(1,ikpt)))
222 [ + + ]: 2357636 : if (abs(kpts(2,ikpt)) > tol10) smallestlen = min(smallestlen, abs(kpts(2,ikpt)))
223 [ + + ]: 2358271 : if (abs(kpts(3,ikpt)) > tol10) smallestlen = min(smallestlen, abs(kpts(3,ikpt)))
224 : : end do
225 : 635 : new%max_linear_density = nint(one/smallestlen)
226 : : else
227 : : ! Get it from input
228 : 478 : new%max_linear_density = max_linear_density
229 : : end if
230 : :
231 : 1113 : new%npoints = nkpt
232 : : new%min_rank = nint(real(new%max_linear_density)*(half+tol8 +&
233 : : real(new%max_linear_density)*(half+tol8 +&
234 : 1113 : real(new%max_linear_density)*(half+tol8))))
235 : :
236 : : new%max_rank = nint(real(new%max_linear_density)*(1+half+tol8 +&
237 : : real(new%max_linear_density)*(1+half+tol8 +&
238 : 1113 : real(new%max_linear_density)*(1+half+tol8))))
239 : :
240 : 1113 : timrev = 2
241 : 1113 : new%time_reversal = .true.
242 [ + + ]: 1113 : if (present(time_reversal)) then
243 [ - + ]: 38 : if (.not. time_reversal) timrev = 1
244 : 38 : new%time_reversal = .false.
245 : : end if
246 : :
247 : : ! Ensure kpt(i)+one is positive, and the smallest difference between kpts should be larger than 1/100 ie ngkpt < 100.
248 : : ! the following fills invrank for the k-points in the list provided (may be only the irred kpts)
249 [ + + ]: 1113 : if (compute_invrank_) then
250 [ - + + - : 1938 : ABI_MALLOC(new%invrank, (new%min_rank:new%max_rank))
+ - + - -
+ - + -
+ ]
251 [ + + ]: 2506059 : new%invrank(:) = -1
252 : :
253 [ + + ]: 2360096 : do ikpt=1,nkpt
254 : 2359450 : irank = new%get_rank(kpts(:,ikpt))
255 [ + - - + ]: 2359450 : if (irank > new%max_rank .or. irank < new%min_rank) then
256 : 0 : write(msg,'(a,2i0)')" rank above max_rank or below min_rank, ikpt, rank ", ikpt, irank
257 : 0 : ABI_ERROR(msg)
258 : : end if
259 : 2360096 : new%invrank(irank) = ikpt
260 : : end do
261 : : end if
262 : :
263 : : ! if symrec is provided, fill invrank with appropriate irred kpt indices
264 : : ! for symmetry completion: kptrank_t%invrank points to the irred k-point
265 : : ! equivalent to the k-point whose rank is provided
266 [ + + ]: 1113 : if (present(symrec)) then
267 [ - + ]: 39 : if(.not. present(nsym)) then
268 : 0 : ABI_ERROR("need both symrec and nsym arguments together")
269 : : end if
270 [ + + ]: 1307 : do ikpt=1,nkpt
271 : : ! itim == 1 for positive, and itim==2 gives Kramers opposite of k-point
272 : : ! favor the former by looping it last
273 [ + + ]: 3843 : do itim = timrev, 1, -1
274 [ + + ]: 123132 : do isym = 1, nsym
275 [ + + + + : 3341184 : symkpt = (-1)**(itim+1) * matmul(symrec(:,:,isym), kpts(:, ikpt))
+ + + + +
+ ]
276 : 119328 : symkptrank = new%get_rank(symkpt(:))
277 : 121864 : new%invrank(symkptrank) = ikpt
278 : : end do
279 : : end do
280 : : end do
281 : : end if
282 : :
283 : 1113 : end subroutine krank_init
284 : : !!***
285 : :
286 : : !----------------------------------------------------------------------
287 : :
288 : : !!****f* m_krank/get_rank
289 : : !! NAME
290 : : !! get_rank
291 : : !!
292 : : !! FUNCTION
293 : : !! Calculate the rank for one kpt.
294 : : !!
295 : : !! INPUTS
296 : : !! kpt = coordinates of kpoints
297 : : !!
298 : : !! OUTPUT
299 : : !! rank = rank of the kpoint
300 : : !!
301 : : !! SOURCE
302 : :
303 : 764648081 : integer function get_rank(krank, kpt) result(rank)
304 : :
305 : : !Arguments ------------------------------------
306 : : !scalars
307 : : class(krank_t), intent(in) :: krank
308 : : !arrays
309 : : real(dp),intent(in) :: kpt(3)
310 : :
311 : : !Local variables-------------------------------
312 : : !scalars
313 : : character(len=500) :: msg
314 : : !arrays
315 : : real(dp) :: redkpt(3)
316 : : ! *************************************************************************
317 : :
318 : : ! wrap to [0, 1[ -> replaced call to wrap2_zero2one inline, to encapsulate this module
319 [ + + ]: 764648081 : if (kpt(1)>zero) then
320 : 329335175 : redkpt(1)=mod((kpt(1)+tol12),one)-tol12
321 : : else
322 : 435312906 : redkpt(1)=-mod(-(kpt(1)-one+tol12),one)+one-tol12
323 : : end if
324 [ + + ]: 764648081 : if(abs(redkpt(1))<tol12)redkpt(1)=zero
325 : :
326 [ + + ]: 764648081 : if (kpt(2)>zero) then
327 : 329193285 : redkpt(2)=mod((kpt(2)+tol12),one)-tol12
328 : : else
329 : 435454796 : redkpt(2)=-mod(-(kpt(2)-one+tol12),one)+one-tol12
330 : : end if
331 [ + + ]: 764648081 : if(abs(redkpt(2))<tol12)redkpt(2)=zero
332 : :
333 [ + + ]: 764648081 : if (kpt(3)>zero) then
334 : 329460858 : redkpt(3)=mod((kpt(3)+tol12),one)-tol12
335 : : else
336 : 435187223 : redkpt(3)=-mod(-(kpt(3)-one+tol12),one)+one-tol12
337 : : end if
338 [ + + ]: 764648081 : if(abs(redkpt(3))<tol12)redkpt(3)=zero
339 : :
340 : : ! rank = int(real(krank%max_linear_density)*(redkpt(3)+half+tol8 + &
341 : : ! real(krank%max_linear_density)*(redkpt(2)+half+tol8 + &
342 : : ! real(krank%max_linear_density)*(redkpt(1)+half+tol8))))
343 : : rank = nint(real(krank%max_linear_density)*(redkpt(1)+half+tol8 + &
344 : : real(krank%max_linear_density)*(redkpt(2)+half+tol8 + &
345 : 764648081 : real(krank%max_linear_density)*(redkpt(3)+half+tol8))))
346 : :
347 [ - + ]: 764648081 : if (rank > krank%max_rank) then
348 : 0 : write(msg,'(2(a,i0))') ' Rank should be <= max_rank: ', krank%max_rank, ' but got: ', rank
349 : 0 : ABI_ERROR(msg)
350 : : end if
351 [ - + ]: 764648081 : if (rank < krank%min_rank) then
352 : : !print *, "redkpt", redkpt
353 : 0 : write(msg,'(2(a,i0))') ' Rank should be >= min_rank ', krank%min_rank, ' but got: ', rank
354 : 0 : ABI_ERROR(msg)
355 : : end if
356 : :
357 : 764648081 : end function get_rank
358 : : !!***
359 : :
360 : : !----------------------------------------------------------------------
361 : :
362 : : !!****f* m_krank/krank_get_index
363 : : !! NAME
364 : : !! krank_get_index
365 : : !!
366 : : !! FUNCTION
367 : : !! Return the index of the k-point `kpt` in the initial set. -1 if not found.
368 : : !!
369 : : !! INPUTS
370 : : !! krank = rank object for the k-grid we are using
371 : : !! kpt = coordinates of kpoints
372 : : !!
373 : : !! OUTPUT
374 : : !!
375 : : !! SOURCE
376 : :
377 : 264468193 : integer function krank_get_index(krank, kpt) result(ikpt)
378 : :
379 : : !Arguments ------------------------------------
380 : : !scalars
381 : : class(krank_t), intent(in) :: krank
382 : : !arrays
383 : : real(dp),intent(in) :: kpt(3)
384 : :
385 : : !Local variables-------------------------------
386 : : integer :: kpt_rank
387 : : ! *************************************************************************
388 : :
389 : 264468193 : kpt_rank = krank%get_rank(kpt)
390 : 264468193 : ikpt = -1
391 [ + - ]: 264468193 : if (kpt_rank < krank%max_rank) ikpt = krank%invrank(kpt_rank)
392 : :
393 : 264468193 : end function krank_get_index
394 : : !!***
395 : :
396 : : !----------------------------------------------------------------------
397 : :
398 : : !!****f* m_krank/krank_copy
399 : : !!
400 : : !! NAME
401 : : !! krank_copy
402 : : !!
403 : : !! FUNCTION
404 : : !! Deep copy of the object
405 : : !!
406 : : !! INPUTS
407 : : !!
408 : : !! OUTPUT
409 : : !! krank = object containing ranking and inverse ranking, to be deallocated
410 : : !!
411 : : !! SOURCE
412 : :
413 : 33 : type(krank_t) function krank_copy(krank_in) result(krank_out)
414 : :
415 : : !Arguments ------------------------------------
416 : : class(krank_t), intent(in) :: krank_in
417 : : ! *********************************************************************
418 : :
419 : 33 : krank_out%max_linear_density = krank_in%max_linear_density
420 : 33 : krank_out%min_rank = krank_in%min_rank
421 : 33 : krank_out%max_rank = krank_in%max_rank
422 : 33 : krank_out%npoints = krank_in%npoints
423 : :
424 [ + - + - : 99 : ABI_MALLOC(krank_out%invrank, (krank_out%min_rank:krank_out%max_rank))
+ - - + -
+ ]
425 [ + - - + : 8069 : krank_out%invrank = krank_in%invrank
- - - - +
+ ]
426 : :
427 : : ! This is why I call it deep copy!
428 [ - + + - : 99 : ABI_MALLOC(krank_out%kpts, (3, size(krank_in%kpts, dim=2)))
+ - + - -
+ - + ]
429 [ + + + + ]: 12697 : krank_out%kpts = krank_in%kpts
430 : 33 : krank_out%kpts_owns_memory = .True.
431 : :
432 [ - + ]: 33 : end function krank_copy
433 : : !!***
434 : :
435 : : !----------------------------------------------------------------------
436 : :
437 : : !!****f* m_krank/krank_free
438 : : !! NAME
439 : : !! krank_free
440 : : !!
441 : : !! FUNCTION
442 : : !! This routine deallocates the arrays in a krank_t structure
443 : : !!
444 : : !! INPUTS
445 : : !! krank = object containing ranking and inverse ranking, to be deallocated
446 : : !!
447 : : !! SOURCE
448 : :
449 : 1209 : subroutine krank_free(krank)
450 : :
451 : : !Arguments ------------------------------------
452 : : class(krank_t), intent(inout) :: krank
453 : : ! *********************************************************************
454 : :
455 [ + + ]: 1209 : ABI_SFREE(krank%invrank)
456 [ + + ]: 1209 : ABI_SFREE(krank%rank2ikpt_)
457 [ + + ]: 1209 : ABI_SFREE(krank%rank2symtime_)
458 : :
459 [ + + ]: 1209 : if (krank%kpts_owns_memory) then
460 [ + - ]: 33 : ABI_SFREE_PTR(krank%kpts)
461 : 33 : krank%kpts_owns_memory = .False.
462 : 33 : krank%kpts => null()
463 : : else
464 : 1176 : krank%kpts => null()
465 : : end if
466 : :
467 : 1209 : end subroutine krank_free
468 : : !!***
469 : :
470 : : !----------------------------------------------------------------------
471 : :
472 : : !!****f* m_krank/krank_print
473 : : !!
474 : : !! NAME
475 : : !! krank_print
476 : : !!
477 : : !! FUNCTION
478 : : !! This routine prints the arrays and dimensions of a krank_t structure
479 : : !!
480 : : !! INPUTS
481 : : !! krank = object containing ranking and inverse ranking
482 : : !! unout = unit for open file to print to
483 : : !!
484 : : !! SOURCE
485 : :
486 : 0 : subroutine krank_print(krank, unout)
487 : :
488 : : !Arguments ------------------------------------
489 : : !scalars
490 : : integer, intent(in) :: unout
491 : : !arrays
492 : : class(krank_t), intent(in) :: krank
493 : : ! *********************************************************************
494 : :
495 : 0 : write(unout, *)
496 : 0 : write(unout, '(a)') ' Dump of the contents of a krank_t structure with k-point rank information'
497 : 0 : write(unout, '(a,i0)') ' max linear density of points in 3 directions: max_linear_density = ', krank%max_linear_density
498 : 0 : write(unout, '(a,i0)') ' maximum rank for any point in grid: max_rank = ', krank%max_rank
499 : 0 : write(unout, '(a,i0)') ' number of points in input grid: npoints = ', krank%npoints
500 : 0 : write(unout, *)
501 : 0 : write(unout, '(a)') ' invrank array = '
502 : 0 : write(unout, '(i0)') krank%invrank(:)
503 : 0 : write(unout, *)
504 : :
505 : 0 : end subroutine krank_print
506 : : !!***
507 : :
508 : : !----------------------------------------------------------------------
509 : :
510 : : !!****f* m_krank/krank_get_mapping
511 : : !! NAME
512 : : !! krank_get_mapping
513 : : !!
514 : : !! FUNCTION
515 : : !! Use symmetries to map input kptn2 to the list of k-points used to generate krank_t.
516 : : !! Similar to listkk but, unlike listkk, this algo does not try to minimize the distance
517 : : !! Mainly used to map two set of k-points associated to the same grid (e.g. BZ --> IBZ, IBZ(q) --> IBZ etc.
518 : : !! Must be faster than listkk for dense meshes
519 : : !! although this routine requires the allocation of temporary array of shape (2, self%min_rank:self%max_rank)
520 : : !! Returns indirect indexing list indkk.
521 : : !!
522 : : !! INPUTS
523 : : !! kptns2(3,nkpt2)=list of final k points
524 : : !! nkpt2=number of final k points
525 : : !! gmet(3,3)=reciprocal space metric (bohr^-2)
526 : : !! nsym=number of symmetry elements
527 : : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
528 : : !! symmat(3,3,nsym)=symmetry operations (symrel or symrec, depending on value of use_symrec)
529 : : !! timrev=1 if the use of time-reversal is allowed; 0 otherwise
530 : : !! [use_symrec]: if present and true, symmat assumed to be symrec, otherwise assumed to be symrel (default)
531 : : !! [qpt]: q-point to be added to kptns2. 0 if not specified.
532 : : !!
533 : : !! OUTPUT
534 : : !! dksqmax=maximal value of the norm**2 of the difference between
535 : : !! a kpt2 vector and the closest k-point found from the kptns1 set, using symmetries.
536 : : !! indkk(6, nkpt2)=
537 : : !! indkk(:,1)=k point index of kptns1
538 : : !! indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
539 : : !! indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
540 : : !! to give kpt1b, that is the closest to kpt2.
541 : : !! indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
542 : : !!
543 : : !! SOURCE
544 : :
545 : 87045 : subroutine krank_get_mapping(self, nkpt2, kptns2, dksqmax, gmet, indkk, nsym, symafm, symmat, timrev, &
546 : : use_symrec, qpt) ! optional
547 : :
548 : : !Arguments ------------------------------------
549 : : !scalars
550 : : class(krank_t),intent(inout) :: self
551 : : integer,intent(in) :: nkpt2, nsym, timrev
552 : : real(dp),intent(out) :: dksqmax
553 : : logical,optional,intent(in) :: use_symrec
554 : : !arrays
555 : : integer,intent(in) :: symafm(nsym), symmat(3,3,nsym)
556 : : integer,intent(out) :: indkk(6, nkpt2)
557 : : real(dp),intent(in) :: gmet(3,3), kptns2(3,nkpt2)
558 : : real(dp),optional,intent(in) :: qpt(3)
559 : :
560 : : !Local variables-------------------------------
561 : : !scalars
562 : : integer :: irank, ikpt1, ikpt2, itimrev, isym, ii
563 : : logical :: my_use_symrec
564 : : !arrays
565 : 174090 : integer :: dkint(3), my_symmat(3, 3, nsym)
566 : : !integer,allocatable :: rank2ikpt(:), rank2symtime(:)
567 : : real(dp) :: kpt1a(3), dk(3), my_qpt(3)
568 : : ! *************************************************************************
569 : :
570 : : !print *,'self%kpts (gstore%ibz)'
571 : : !print *, self%kpts
572 : : !print *,'kptns2 (gstore%bz)'
573 : : !print *, kptns2
574 : : !print *
575 : :
576 [ + - ]: 87045 : my_use_symrec = .False.; if (present(use_symrec)) my_use_symrec = use_symrec
577 [ + - ]: 87045 : my_qpt = zero; if (present(qpt)) my_qpt = qpt
578 : :
579 [ + + ]: 87045 : if (my_use_symrec) then
580 : : ! Symrec k
581 [ + + + + : 491355 : my_symmat = symmat
+ + ]
582 : : else
583 : : ! Symrel^T k
584 [ + + ]: 4204139 : do isym=1,nsym
585 [ + + + + ]: 53619023 : my_symmat(:,:,isym) = transpose(symmat(:,:,isym))
586 : : !print *, 'isym =', isym
587 : : !print *, my_symmat(1, :, isym)
588 : : !print *, my_symmat(2, :, isym)
589 : : !print *, my_symmat(3, :, isym)
590 : : !print *
591 : : end do
592 : : end if
593 : :
594 [ + + - + : 87961 : ABI_MALLOC_IFNOT(self%rank2symtime_, (self%min_rank:self%max_rank))
+ - + - +
- - + - +
- + ]
595 [ + + - + : 87961 : ABI_MALLOC_IFNOT(self%rank2ikpt_, (self%min_rank:self%max_rank))
+ - + - +
- - + - +
- + ]
596 [ + + ]: 141473081 : self%rank2ikpt_ = -1 !; self%rank2symtime_ = -1
597 : :
598 [ + + ]: 5219396 : do ikpt1=1,self%npoints
599 : :
600 [ + + ]: 15484098 : do itimrev=0,timrev
601 [ + + ]: 505992845 : do isym=1,nsym
602 : : ! Do not use magnetic symmetries.
603 [ + - ]: 490595792 : if (symafm(isym) == -1) cycle
604 : :
605 [ + + + + : 13736682176 : kpt1a = (1 - 2*itimrev) * matmul(my_symmat(:, :, isym), self%kpts(:, ikpt1))
+ + + + +
- + + ]
606 : 490595792 : irank = self%get_rank(kpt1a)
607 [ + + ]: 500860494 : if (self%rank2ikpt_(irank) == -1) then
608 : 130023831 : self%rank2ikpt_(irank) = ikpt1
609 : 130023831 : self%rank2symtime_(irank) = isym + itimrev * nsym
610 : : end if
611 : : end do
612 : : end do
613 : :
614 : : end do
615 : :
616 : 87045 : dksqmax = zero
617 [ + + ]: 5453834 : do ikpt2=1,nkpt2
618 [ + + ]: 21467156 : irank = self%get_rank(kptns2(:, ikpt2) + my_qpt)
619 : 5366789 : ikpt1 = self%rank2ikpt_(irank)
620 : 5366789 : ii = self%rank2symtime_(irank)
621 : 5366789 : isym = 1 + mod(ii - 1, nsym)
622 : 5366789 : itimrev = (ii - 1) / nsym
623 [ + + + + : 150270092 : kpt1a = (1 - 2 * itimrev) * matmul(my_symmat(:, :, isym), self%kpts(:, ikpt1))
+ + + + +
- + + ]
624 [ + + ]: 21467156 : dk(:) = kptns2(:,ikpt2) + my_qpt - kpt1a(:)
625 [ + + ]: 21467156 : dkint(:) = nint(dk(:) + tol12)
626 : :
627 : : !print *, 'BZ point', ikpt2 !'kptns2(:, ', ikpt2, ')'
628 : : !print *, kptns2(:, ikpt2)
629 : : !print *, 'isym =', isym
630 : : !print *, my_symmat(1, :, isym)
631 : : !print *, my_symmat(2, :, isym)
632 : : !print *, my_symmat(3, :, isym)
633 : : !print *, 'time_reversal =', itimrev
634 : : !print *, 'IBZ point', ikpt1 !'self%kpts(:, ', ikpt1, ')'
635 : : !!print *, self%kpts(:, ikpt1)
636 : : !!print *, 'kpt1a'
637 : : !!print *, kpt1a
638 : : !print *, 'difference % lattice' !'dk'
639 : : !print *, dk
640 : : !print *
641 : :
642 : 5366789 : indkk(1, ikpt2) = ikpt1
643 : 5366789 : indkk(2, ikpt2) = isym
644 [ + + ]: 21467156 : indkk(3:5, ikpt2) = dkint(:)
645 : 5366789 : indkk(6, ikpt2) = itimrev
646 : :
647 : : ! Compute norm of the difference vector.
648 [ + + ]: 21467156 : dk(:) = dk(:) - dkint(:)
649 : :
650 : : dksqmax = max(dksqmax, &
651 : : gmet(1,1)*dk(1)**2 + gmet(2,2)*dk(2)**2 + gmet(3,3)*dk(3)**2 + &
652 : 5453834 : two * (gmet(2,1)*dk(2)*dk(1) + gmet(3,2)*dk(3)*dk(2)+gmet(3,1)*dk(3)*dk(1)))
653 : :
654 : : !if (dksqmax > tol8) then
655 : : ! print *, "kbase:", self%kpts(:, ikpt1)
656 : : ! print *, "k", kptns2(:, ikpt2)
657 : : ! print *, "k + q", kptns2(:, ikpt2) + my_qpt
658 : : ! print *, "dk", dk
659 : : ! print *, "dksqmax", dksqmax
660 : : ! print *
661 : : !end if
662 : : end do
663 : :
664 : : !ABI_FREE(self%rank2ikpt_)
665 : : !ABI_FREE(self%rank2symtime_)
666 : :
667 : 87045 : end subroutine krank_get_mapping
668 : : !!***
669 : :
670 : : !----------------------------------------------------------------------
671 : :
672 : : !!****f* m_krank/get_ibz2bz
673 : : !! NAME
674 : : !! get_ibz2bz
675 : : !!
676 : : !! FUNCTION
677 : : !! Return array with the index of the IBZ wave vectors in the BZ.
678 : : !!
679 : : !! INPUTS
680 : : !!
681 : : !! OUTPUT
682 : : !!
683 : : !! SOURCE
684 : :
685 : 16 : subroutine get_ibz2bz(nibz, nbz, bz2ibz, ibz2bz, err_msg, ierr)
686 : :
687 : : !Arguments ------------------------------------
688 : : !scalars
689 : : integer,intent(in) :: nibz, nbz
690 : : integer,intent(in) :: bz2ibz(6, nbz)
691 : : integer,intent(out) :: ierr
692 : : character(len=*),intent(out) :: err_msg
693 : : !arrays
694 : : integer,allocatable,intent(out) :: ibz2bz(:)
695 : :
696 : : !Local variables-------------------------------
697 : : !scalars
698 : : integer :: iq_bz, iq_ibz, isym_q, trev_q, cnt, g0_q(3)
699 : : logical :: isirr_q
700 : : !----------------------------------------------------------------------
701 : :
702 [ + - + - : 48 : ABI_MALLOC(ibz2bz, (nibz))
+ - - + -
+ - + ]
703 : :
704 : 16 : cnt = 0
705 [ + + ]: 5130 : do iq_bz=1,nbz
706 : 5114 : iq_ibz = bz2ibz(1, iq_bz); isym_q = bz2ibz(2, iq_bz)
707 [ + + ]: 20456 : trev_q = bz2ibz(6, iq_bz); g0_q = bz2ibz(3:5,iq_bz)
708 [ + + + + : 17492 : isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0))
+ + + - ]
709 : 16 : if (isirr_q) then
710 : 266 : cnt = cnt + 1
711 : 266 : ibz2bz(iq_ibz) = iq_bz
712 : : end if
713 : : end do
714 : :
715 [ - + ]: 16 : ierr = merge(0, 1, cnt == nibz)
716 [ + - ]: 16 : err_msg = ""
717 [ - + ]: 16 : if (ierr /= 0) then
718 [ # # # # ]: 0 : err_msg = sjoin("The number of points in the IBZ computed from symmetry table is: ", itoa(cnt), " while it should be: ", itoa(nibz))
719 : : end if
720 : :
721 [ - + ]: 16 : end subroutine get_ibz2bz
722 : : !!***
723 : :
724 : : !----------------------------------------------------------------------
725 : :
726 : : !!****f* m_krank/star_from_ibz_idx
727 : : !! NAME
728 : : !! star_from_ibz_idx
729 : : !!
730 : : !! FUNCTION
731 : : !! Return array with the indices of the star of the ik_ibz wavevector in the IBZ.
732 : : !! Return number of points in the star, allocate array with the indices of the
733 : : !! star points in the BZ.
734 : : !!
735 : : !! INPUTS
736 : : !!
737 : : !! OUTPUT
738 : : !!
739 : : !! SOURCE
740 : :
741 : 8 : subroutine star_from_ibz_idx(ik_ibz, nkbz, bz2ibz, nk_in_star, kstar_bz_inds)
742 : :
743 : : !Arguments ------------------------------------
744 : : integer,intent(in) :: ik_ibz, nkbz
745 : : integer,intent(out) :: nk_in_star
746 : : integer,intent(in) :: bz2ibz(6, nkbz)
747 : : integer,allocatable,intent(out) :: kstar_bz_inds(:)
748 : :
749 : : !Local variables-------------------------------
750 : : !scalars
751 : : integer :: iq_bz
752 : :
753 : : !----------------------------------------------------------------------
754 : :
755 [ + + + + ]: 17584 : nk_in_star = count(bz2ibz(1, :) == ik_ibz)
756 [ - + + - : 24 : ABI_MALLOC(kstar_bz_inds, (nk_in_star))
+ - + - -
+ - + -
+ ]
757 : :
758 : 8 : nk_in_star = 0
759 [ + + ]: 17584 : do iq_bz=1,nkbz
760 [ + + ]: 17576 : if (bz2ibz(1, iq_bz) /= ik_ibz) cycle
761 : 89 : nk_in_star = nk_in_star + 1
762 : 17584 : kstar_bz_inds(nk_in_star) = iq_bz
763 : : end do
764 : :
765 : 8 : end subroutine star_from_ibz_idx
766 : : !!***
767 : :
768 : : !----------------------------------------------------------------------
769 : :
770 [ - + - - : 4773 : end module m_krank
+ + - + +
- - + - +
+ - - + ]
771 : : !!***
|
