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 87045 : my_use_symrec = .False.; if (present(use_symrec)) my_use_symrec = use_symrec
571 87045 : my_qpt = zero; if (present(qpt)) my_qpt = qpt
572 :
573 87045 : if (my_use_symrec) then
574 : ! Symrec k
575 491355 : my_symmat = symmat
576 : else
577 : ! Symrel^T k
578 4204139 : do isym=1,nsym
579 53619023 : my_symmat(:,:,isym) = transpose(symmat(:,:,isym))
580 : end do
581 : end if
582 :
583 87961 : ABI_MALLOC_IFNOT(self%rank2symtime_, (self%min_rank:self%max_rank))
584 87961 : ABI_MALLOC_IFNOT(self%rank2ikpt_, (self%min_rank:self%max_rank))
585 141473081 : self%rank2ikpt_ = -1 !; self%rank2symtime_ = -1
586 :
587 5219396 : do ikpt1=1,self%npoints
588 :
589 15484098 : do itimrev=0,timrev
590 505992845 : do isym=1,nsym
591 : ! Do not use magnetic symmetries.
592 490595792 : if (symafm(isym) == -1) cycle
593 :
594 13736682176 : kpt1a = (1 - 2*itimrev) * matmul(my_symmat(:, :, isym), self%kpts(:, ikpt1))
595 490595792 : irank = self%get_rank(kpt1a)
596 500860494 : if (self%rank2ikpt_(irank) == -1) then
597 130023831 : self%rank2ikpt_(irank) = ikpt1
598 130023831 : self%rank2symtime_(irank) = isym + itimrev * nsym
599 : end if
600 : end do
601 : end do
602 :
603 : end do
604 :
605 87045 : dksqmax = zero
606 5453834 : do ikpt2=1,nkpt2
607 21467156 : irank = self%get_rank(kptns2(:, ikpt2) + my_qpt)
608 5366789 : ikpt1 = self%rank2ikpt_(irank)
609 5366789 : ii = self%rank2symtime_(irank)
610 5366789 : isym = 1 + mod(ii - 1, nsym)
611 5366789 : itimrev = (ii - 1) / nsym
612 150270092 : kpt1a = (1 - 2 * itimrev) * matmul(my_symmat(:, :, isym), self%kpts(:, ikpt1))
613 21467156 : dk(:) = kptns2(:,ikpt2) + my_qpt - kpt1a(:)
614 21467156 : dkint(:) = nint(dk(:) + tol12)
615 :
616 5366789 : indkk(1, ikpt2) = ikpt1
617 5366789 : indkk(2, ikpt2) = isym
618 21467156 : indkk(3:5, ikpt2) = dkint(:)
619 5366789 : indkk(6, ikpt2) = itimrev
620 :
621 : ! Compute norm of the difference vector.
622 21467156 : dk(:) = dk(:) - dkint(:)
623 :
624 : dksqmax = max(dksqmax, &
625 : gmet(1,1)*dk(1)**2 + gmet(2,2)*dk(2)**2 + gmet(3,3)*dk(3)**2 + &
626 5453834 : two * (gmet(2,1)*dk(2)*dk(1) + gmet(3,2)*dk(3)*dk(2)+gmet(3,1)*dk(3)*dk(1)))
627 :
628 : !if (dksqmax > tol8) then
629 : ! print *, "kbase:", self%kpts(:, ikpt1)
630 : ! print *, "k", kptns2(:, ikpt2)
631 : ! print *, "k + q", kptns2(:, ikpt2) + my_qpt
632 : ! print *, "dk", dk, "dksqmax", dksqmax
633 : !end if
634 : end do
635 :
636 : !ABI_FREE(self%rank2ikpt_)
637 : !ABI_FREE(self%rank2symtime_)
638 :
639 87045 : end subroutine krank_get_mapping
640 : !!***
641 :
642 : !----------------------------------------------------------------------
643 :
644 : !!****f* m_krank/get_ibz2bz
645 : !! NAME
646 : !! get_ibz2bz
647 : !!
648 : !! FUNCTION
649 : !! Return array with the index of the IBZ wave vectors in the BZ.
650 : !!
651 : !! INPUTS
652 : !!
653 : !! OUTPUT
654 : !!
655 : !! SOURCE
656 :
657 16 : subroutine get_ibz2bz(nibz, nbz, bz2ibz, ibz2bz, err_msg, ierr)
658 :
659 : !Arguments ------------------------------------
660 : !scalars
661 : integer,intent(in) :: nibz, nbz
662 : integer,intent(in) :: bz2ibz(6, nbz)
663 : integer,intent(out) :: ierr
664 : character(len=*),intent(out) :: err_msg
665 : !arrays
666 : integer,allocatable,intent(out) :: ibz2bz(:)
667 :
668 : !Local variables-------------------------------
669 : !scalars
670 : integer :: iq_bz, iq_ibz, isym_q, trev_q, cnt, g0_q(3)
671 : logical :: isirr_q
672 : !----------------------------------------------------------------------
673 :
674 48 : ABI_MALLOC(ibz2bz, (nibz))
675 :
676 16 : cnt = 0
677 5130 : do iq_bz=1,nbz
678 5114 : iq_ibz = bz2ibz(1, iq_bz); isym_q = bz2ibz(2, iq_bz)
679 20456 : trev_q = bz2ibz(6, iq_bz); g0_q = bz2ibz(3:5,iq_bz)
680 17492 : isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0))
681 16 : if (isirr_q) then
682 266 : cnt = cnt + 1
683 266 : ibz2bz(iq_ibz) = iq_bz
684 : end if
685 : end do
686 :
687 16 : ierr = merge(0, 1, cnt == nibz)
688 16 : err_msg = ""
689 16 : if (ierr /= 0) then
690 0 : err_msg = sjoin("The number of points in the IBZ computed from symmetry table is: ", itoa(cnt), " while it should be: ", itoa(nibz))
691 : end if
692 :
693 16 : end subroutine get_ibz2bz
694 : !!***
695 :
696 : !----------------------------------------------------------------------
697 :
698 : !!****f* m_krank/star_from_ibz_idx
699 : !! NAME
700 : !! star_from_ibz_idx
701 : !!
702 : !! FUNCTION
703 : !! Return array with the indices of the star of the ik_ibz wavevector in the IBZ.
704 : !! Return number of points in the star, allocate array with the indices of the
705 : !! star points in the BZ.
706 : !!
707 : !! INPUTS
708 : !!
709 : !! OUTPUT
710 : !!
711 : !! SOURCE
712 :
713 8 : subroutine star_from_ibz_idx(ik_ibz, nkbz, bz2ibz, nk_in_star, kstar_bz_inds)
714 :
715 : !Arguments ------------------------------------
716 : integer,intent(in) :: ik_ibz, nkbz
717 : integer,intent(out) :: nk_in_star
718 : integer,intent(in) :: bz2ibz(6, nkbz)
719 : integer,allocatable,intent(out) :: kstar_bz_inds(:)
720 :
721 : !Local variables-------------------------------
722 : !scalars
723 : integer :: iq_bz
724 :
725 : !----------------------------------------------------------------------
726 :
727 17584 : nk_in_star = count(bz2ibz(1, :) == ik_ibz)
728 24 : ABI_MALLOC(kstar_bz_inds, (nk_in_star))
729 :
730 8 : nk_in_star = 0
731 17584 : do iq_bz=1,nkbz
732 17576 : if (bz2ibz(1, iq_bz) /= ik_ibz) cycle
733 89 : nk_in_star = nk_in_star + 1
734 17584 : kstar_bz_inds(nk_in_star) = iq_bz
735 : end do
736 :
737 8 : end subroutine star_from_ibz_idx
738 : !!***
739 :
740 : !----------------------------------------------------------------------
741 :
742 4773 : end module m_krank
743 : !!***
|
