Line data Source code
1 : !!****m* ABINIT/m_kpts
2 : !! NAME
3 : !! m_kpts
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2008-2026 ABINIT group (XG, MG, MJV, DRH, DCA, JCC, MM)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_kpts
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 : use m_crystal
27 : use m_sort
28 : use m_krank
29 : use m_htetra
30 : use m_xmpi
31 :
32 : use m_time, only : timab, cwtime, cwtime_report
33 : use m_copy, only : alloc_copy
34 : use m_numeric_tools, only : wrap2_zero_one, interpol3d_0d
35 : use m_matrix, only : mati3inv, mati3det, matr3inv
36 : use m_symtk, only : smallprim
37 : use m_fstrings, only : sjoin, itoa, ftoa, ltoa, ktoa
38 : use m_numeric_tools, only : wrap2_pmhalf
39 : use m_geometry, only : metric
40 : use m_symkpt, only : symkpt, symkpt_new
41 :
42 : implicit none
43 :
44 : private
45 :
46 : public :: kpts_timrev_from_kptopt ! Returns the value of timrev from kptopt
47 : public :: kpts_ibz_from_kptrlatt ! Determines the IBZ, the weights and the BZ from kptrlatt
48 : public :: tetra_from_kptrlatt ! Create an instance of htetra_t from kptrlatt and shiftk
49 : public :: symkchk ! Checks that the set of k points has the full space group symmetry, modulo time reversal if appropriate.
50 : public :: kpts_sort ! order list of k-points according to the norm.
51 : public :: kpts_pack_in_stars ! Pack k-points in stars.
52 : public :: kpts_map ! Compute symmetry table.
53 : public :: kpts_map_print ! Print the symmetry table bz2ibz to a list of units with header.
54 : public :: listkk ! Find correspondence between two set of k-points.
55 : public :: getkgrid ! Compute the grid of k points in the irreducible Brillouin zone.
56 : public :: smpbz ! Generate a set of special k (or q) points which samples in a homogeneous way the BZ
57 : public :: testkgrid ! Test different grids of k points.
58 : public :: kptrlatt_from_ngkpt ! Insert ngkpt in kptrlatt matrix
59 :
60 : !FIXME: Deprecated
61 : public :: get_full_kgrid ! Create full grid of kpoints and find equivalent irred ones.
62 : public :: mknormpath
63 : private :: get_kpt_fullbz ! Create full grid of kpoints from kptrlatt and shiftk
64 : !!***
65 :
66 : #if 0
67 : type, public :: kinds_t
68 : integer :: ibz_idx
69 : integer :: isym
70 : integer :: trev
71 : integer :: g0(3)
72 : logical :: is_irred
73 : contains
74 : end type kinds_t
75 : !!***
76 :
77 : type, public :: kmap_t
78 : integer :: nkpt = -1
79 : integer :: kptopt = -1
80 : integer :: ierr = -1
81 : !character(len=500) :: err_msg
82 : character(len=20) :: mode
83 : real(dp) :: dksqmax
84 : real(dp) :: qpt(3) = zero
85 : type(kinds_t),allocatable :: inds(:)
86 : contains
87 : procedure :: free => kmap_free
88 : end type kmap_t
89 : !!***
90 : #endif
91 :
92 : !----------------------------------------------------------------------
93 :
94 : !!****t* m_kpts/bzlint_t
95 : !! NAME
96 : !! bzlint_t
97 : !!
98 : !! FUNCTION
99 : !! Linear interpolator for functions defined in the BZ.
100 : !!
101 : !! To interpolate initial data defined on a ngkpt(3) k-mesh.
102 : !! if the values of shape (ndat, nkpt) are known on nkpt k-points `kpts(3,nkpt)
103 : !! belonging to the ngkpt mesh, use the following calls:
104 : !!
105 : !! Example:
106 : !! type(bzlint_t) :: bzlint
107 : !! call bzlint%init(ngkpt, ndat, nkpt, kpts, values)
108 : !!
109 : !! ! Now we can (linearly) interpolate at arbitrary kpoints
110 : !! allocate(results, (ndat)
111 : !! do ik=1,nk_interp(interp_kpt, results)
112 : !! call bzlint%free(results)
113 : !! end do
114 : !!
115 : !! call bzlint%free() ! Free memory
116 : !!
117 : !! To handle complex data, use real(dp) pointers associated to complex(dp) arrays as in:
118 : !!
119 : !! use, intrinsic :: iso_c_binding
120 : !! complex(dp),allocatable,target :: cvalues(:)
121 : !! real(dp), ABI_CONTIGUOUS pointer :: rpt_d2(:,:)
122 : !!
123 : !! allocate(cvalues(ndat, nkpt))
124 : !! ! fill cvalues...
125 : !!
126 : !! call c_f_pointer(c_loc(cvalues), rpt_d2, [2*ndat, nkpt])
127 : !! call bzlint%init(ngkpt, 2*ndat, nkpt, kpts, rpt_d2)
128 : !!
129 : !! SOURCE
130 :
131 : type, public :: bzlint_t
132 : integer :: nx, ny, nz, ndat
133 : integer :: ngkpt(3)
134 : real(dp),allocatable :: vals_grid(:,:,:,:)
135 : contains
136 : procedure :: init => bzlint_init
137 : procedure :: interp => bzlint_interp
138 : procedure :: free => bzlint_free
139 : end type bzlint_t
140 : !!***
141 :
142 : !----------------------------------------------------------------------
143 :
144 : contains !============================================================
145 : !!***
146 :
147 : !!****f* m_kpts/kpts_timrev_from_kptopt
148 : !! NAME
149 : !! kpts_timrev_from_kptopt
150 : !!
151 : !! FUNCTION
152 : !! Returns the value of timrev from kptopt
153 : !! 1 if the use of time-reversal is allowed; 0 otherwise
154 : !!
155 : !! INPUTS
156 : !! kptopt=option for the generation of k points
157 : !! (defines whether spatial symmetries and/or time-reversal can be used)
158 : !!
159 : !! SOURCE
160 :
161 91099 : integer pure function kpts_timrev_from_kptopt(kptopt) result(timrev)
162 :
163 : !Arguments ------------------------------------
164 : integer,intent(in) :: kptopt
165 : ! *********************************************************************
166 :
167 91099 : timrev = 1; if (any(kptopt == [3, 4])) timrev = 0
168 :
169 4054 : end function kpts_timrev_from_kptopt
170 : !!***
171 :
172 : !!****f* m_kpts/kpts_ibz_from_kptrlatt
173 : !! NAME
174 : !! kpts_ibz_from_kptrlatt
175 : !!
176 : !! FUNCTION
177 : !! Determines the irreducible wedge, the corresponding weights and the list
178 : !! of k-points in the Brillouin Zone starting from kptrlatt and the set shifts.
179 : !!
180 : !! INPUTS
181 : !! cryst<crystal_t> = crystalline structure with info on symmetries and time-reversal.
182 : !! kptopt=option for the generation of k points (defines whether spatial symmetries and/or time-reversal can be used)
183 : !! kptrlatt(3,3)=integer coordinates of the primitive vectors of the
184 : !! lattice reciprocal to the k point lattice to be generated here
185 : !! If diagonal, the three values are the Monkhorst-Pack usual values, in case of simple cubic.
186 : !! nshiftk= number of shift vectors in the repeated cell
187 : !! shiftk(3,nshiftk) = vectors that will be used to determine the shifts from (0. 0. 0.).
188 : !!
189 : !! OUTPUT
190 : !! nkibz,nkbz = Number of points in IBZ and BZ, respectively.
191 : !! The following arrays are allocated and returned by the routine:
192 : !! kibz(3,nkibz) = k-points in the IBZ.
193 : !! wtk(nkibz) = weights of the k-points in the IBZ (normalized to one).
194 : !! kbz(3,nkbz) = k-points in the BZ.
195 : !! [new_kptrlatt] = New value of kptrlatt returned by getkgrid
196 : !! [new_shiftk(3,new_nshiftk)] = New set of shifts returned by getkgrid
197 : !! [bz2ibz(6,nkbz)]=Mapping BZ --> IBZ
198 : !!
199 : !! SOURCE
200 :
201 529 : subroutine kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, & ! in
202 : nkibz, kibz, wtk, nkbz, kbz, & ! out
203 : new_kptrlatt, new_shiftk, bz2ibz) ! out Optional
204 :
205 : !Arguments ------------------------------------
206 : !scalars
207 : integer,intent(in) :: nshiftk,kptopt
208 : integer,intent(out) :: nkibz,nkbz
209 : type(crystal_t),intent(in) :: cryst
210 : !arrays
211 : integer,intent(in) :: kptrlatt(3,3)
212 : integer,optional,allocatable,intent(out) :: bz2ibz(:,:)
213 : integer,optional,intent(out) :: new_kptrlatt(3,3)
214 : real(dp),intent(in) :: shiftk(3,nshiftk)
215 : real(dp),allocatable,intent(out) :: wtk(:),kibz(:,:),kbz(:,:)
216 : real(dp),optional,allocatable,intent(out) :: new_shiftk(:,:)
217 :
218 : !Local variables-------------------------------
219 : !scalars
220 : integer,parameter :: iout0 = 0, chksymbreak0 = 0, iscf2 = 2
221 : integer :: my_nshiftk
222 : real(dp) :: kptrlen
223 : !arrays
224 : integer,parameter :: vacuum0(3) = [0, 0, 0]
225 : integer :: my_kptrlatt(3,3)
226 529 : integer,allocatable :: indkpt(:),bz2ibz_smap(:,:)
227 : real(dp) :: my_shiftk(3,MAX_NSHIFTK)
228 : ! *********************************************************************
229 :
230 : ! Copy kptrlatt and shifts because getkgrid can change them
231 : ! Be careful as getkgrid expects shiftk(3,MAX_NSHIFTK).
232 529 : ABI_CHECK_IRANGE(nshiftk, 1, MAX_NSHIFTK, "Invalid value of nshiftk")
233 2645 : my_nshiftk = nshiftk; my_shiftk = zero; my_shiftk(:,1:nshiftk) = shiftk
234 529 : my_kptrlatt = kptrlatt
235 :
236 : call getkgrid_low(chksymbreak0,iout0,iscf2,kibz,kptopt,my_kptrlatt,kptrlen,&
237 : cryst%nsym,-1,nkibz,my_nshiftk,cryst%nsym,cryst%rprimd,my_shiftk,cryst%symafm,&
238 529 : cryst%symrel,vacuum0,wtk,indkpt,bz2ibz_smap,fullbz=kbz)
239 :
240 529 : if (present(bz2ibz)) then
241 331 : ABI_MOVE_ALLOC(bz2ibz_smap, bz2ibz)
242 : else
243 198 : ABI_SFREE(bz2ibz_smap)
244 : endif
245 529 : ABI_SFREE(indkpt)
246 :
247 529 : nkbz = size(kbz, dim=2)
248 :
249 : ! Optionally, return new shifts and new_kptrlatt
250 529 : if (present(new_shiftk)) then
251 822 : ABI_MALLOC(new_shiftk, (3, my_nshiftk))
252 1644 : new_shiftk = my_shiftk(:, 1:my_nshiftk)
253 : end if
254 529 : if (present(new_kptrlatt)) new_kptrlatt = my_kptrlatt
255 :
256 : DBG_CHECK(abs(sum(wtk) - one) < tol10, "sum(wtk) != one")
257 :
258 529 : end subroutine kpts_ibz_from_kptrlatt
259 : !!***
260 :
261 : !----------------------------------------------------------------------
262 :
263 : !!****f* m_kpts/tetra_from_kptrlatt
264 : !! NAME
265 : !! tetra_from_kptrlatt
266 : !!
267 : !! FUNCTION
268 : !! Helper function to to create an instance and htetra from kptrlatt and shiftk
269 : !!
270 : !! INPUTS
271 : !! cryst<cryst_t>=Crystalline structure.
272 : !! kptopt=Option for the k-point generation.
273 : !! kptrlatt(3,3)=k-point lattice specification
274 : !! nshiftk= number of shift vectors.
275 : !! shiftk(3,nshiftk)=shift vectors for k point generation
276 : !! nkibz=Number of points in the IBZ
277 : !! kibz(3,nkibz)=Reduced coordinates of the k-points in the IBZ.
278 : !! comm= MPI communicator
279 : !!
280 : !! OUTPUT
281 : !! tetra<htetra_t>=Tetrahedron object, fully initialized if ierr == 0.
282 : !! msg=Error message if ierr /= 0
283 : !! ierr=Exit status
284 : !!
285 : !! SOURCE
286 :
287 81 : type(htetra_t) function tetra_from_kptrlatt( &
288 81 : cryst, kptopt, kptrlatt, nshiftk, shiftk, nkibz, kibz, comm, msg, ierr) result (htetra)
289 :
290 : !Arguments ------------------------------------
291 : !scalars
292 : integer,intent(in) :: kptopt,nshiftk,nkibz,comm
293 : integer,intent(out) :: ierr
294 : character(len=*),intent(out) :: msg
295 : type(crystal_t),intent(in) :: cryst
296 : !arrays
297 : integer,intent(in) :: kptrlatt(3,3)
298 : real(dp),intent(in) :: shiftk(3,nshiftk),kibz(3,nkibz)
299 :
300 : !Local variables-------------------------------
301 : !scalars
302 : integer :: nkfull,my_nkibz,new_nshiftk
303 : character(len=80) :: errorstring
304 : !arrays
305 : integer :: new_kptrlatt(3,3)
306 81 : integer,allocatable :: indkk(:)
307 81 : integer,allocatable :: bz2ibz(:,:)
308 : real(dp) :: rlatt(3,3),klatt(3,3)
309 81 : real(dp),allocatable :: kfull(:,:),my_kibz(:,:),my_wtk(:),new_shiftk(:,:)
310 : ! *************************************************************************
311 :
312 81 : ierr = 0
313 :
314 : ! Refuse only 1 kpoint: the algorithms are no longer valid. DOH!
315 81 : if (nkibz == 1) then
316 0 : msg = 'You need at least 2 kpoints to use the tetrahedron method.'
317 0 : ierr = 1; goto 10
318 : end if
319 81 : if (all(kptrlatt == 0)) then
320 0 : msg = 'Cannot generate tetrahedron because input kptrlatt == 0.'
321 0 : ierr = 1; goto 10
322 : end if
323 81 : if (kptopt <= 0) then
324 0 : msg = sjoin("Cannot generate tetrahedron because input kptopt:", itoa(kptopt))
325 0 : ierr = 1; goto 10
326 : end if
327 :
328 : call kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, &
329 81 : my_nkibz, my_kibz, my_wtk, nkfull, kfull, new_kptrlatt=new_kptrlatt, new_shiftk=new_shiftk, bz2ibz=bz2ibz)
330 :
331 81 : ABI_FREE(my_wtk)
332 81 : new_nshiftk = size(new_shiftk, dim=2)
333 :
334 81 : if (my_nkibz /= nkibz .or. all(my_kibz /= kibz) ) then
335 0 : msg = sjoin("Input nkibz:", itoa(nkibz), "does not agree with computed value:", itoa(my_nkibz))
336 0 : ierr = 1; goto 10
337 : end if
338 :
339 : ! Do not support new_nshiftk > 1: lattice must be decomposed into boxes
340 : ! and this is not always possible (I think) with bizarre shifts
341 : ! normally at this point we have incorporated everything into
342 : ! new_kptrlatt, and only 1 shift is needed (in particular for MP grids).
343 81 : if (new_nshiftk > 1) then
344 : write(msg, "(9a)") &
345 0 : 'Cannot create tetrahedron object...',ch10, &
346 0 : 'Only simple lattices are supported. Action: use nshiftk=1.',ch10, &
347 0 : 'new_shiftk: ', trim(ltoa(reshape(new_shiftk, [3*new_nshiftk]))),ch10, &
348 0 : 'new_kptrlatt: ', trim(ltoa(reshape(new_kptrlatt, [9])))
349 0 : ierr = 2; goto 10
350 : end if
351 :
352 1053 : rlatt = new_kptrlatt; call matr3inv(rlatt, klatt)
353 :
354 243 : ABI_MALLOC(indkk, (nkfull))
355 164424 : indkk(:) = bz2ibz(1, :)
356 81 : ABI_SFREE(bz2ibz)
357 :
358 81 : call htetra%init(indkk, cryst%gprimd, klatt, kfull, nkfull, my_kibz, my_nkibz, ierr, errorstring, comm)
359 81 : if (ierr /= 0) msg = errorstring
360 :
361 : 10 continue
362 :
363 81 : ABI_SFREE(my_kibz)
364 81 : ABI_SFREE(indkk)
365 81 : ABI_SFREE(kfull)
366 81 : ABI_SFREE(new_shiftk)
367 :
368 162 : end function tetra_from_kptrlatt
369 : !!***
370 :
371 : !!****f* m_kpts/symkchk
372 : !! NAME
373 : !! symkchk
374 : !!
375 : !! FUNCTION
376 : !! Checks that the set of k points chosen for a response function
377 : !! calculation has the full space group symmetry, modulo time reversal if appropriate.
378 : !! Returns ierr/=0 with error message if not satisfied
379 : !! Currently used only when strain perturbation is treated. Based on symkpt.
380 : !!
381 : !! INPUTS
382 : !! kptns(3,nkpt)= k vectors in reciprocal space
383 : !! nkpt = number of k-points whose weights are wtk
384 : !! nsym=number of space group symmetries
385 : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
386 : !! timrev: if 1, the time reversal operation has to be taken into account
387 : !! if 0, no time reversal symmetry.
388 : !!
389 : !! OUTPUT
390 : !! msg=Error message if ierr /= 0
391 : !!
392 : !! TODO
393 : !! This version should scale badly with the number of k-points. Replace loops with listkk
394 : !!
395 : !! SOURCE
396 :
397 37 : integer function symkchk(kptns,nkpt,nsym,symrec,timrev,errmsg) result(ierr)
398 :
399 : !Arguments -------------------------------
400 : !scalars
401 : integer,intent(in) :: nkpt,nsym,timrev
402 : character(len=*),intent(out) :: errmsg
403 : !arrays
404 : integer,intent(in) :: symrec(3,3,nsym)
405 : real(dp),intent(in) :: kptns(3,nkpt)
406 :
407 : !Local variables -------------------------
408 : !scalars
409 : integer :: identi,ii,ikpt,ikpt2,imatch,isym,jj,tident
410 : real(dp) :: difk,reduce
411 : character(len=500) :: msg
412 : !arrays
413 : real(dp) :: ksym(3)
414 : ! *********************************************************************
415 :
416 37 : ierr = 0
417 :
418 37 : if(timrev/=1 .and. timrev/=0)then
419 : write(errmsg, '(3a,i0,a)' )&
420 0 : 'timrev should be 0 or 1, while',ch10,&
421 0 : 'it is equal to ',timrev,'.'
422 0 : ierr = 1; return
423 : end if
424 :
425 37 : if(nsym/=1)then
426 : ! Find the identity symmetry operation
427 36 : do isym=1,nsym
428 : tident=1
429 144 : do jj=1,3
430 108 : if(symrec(jj,jj,isym)/=1)tident=0
431 468 : do ii=1,3
432 756 : if( ii/=jj .and.symrec(ii,jj,isym)/=0)tident=0
433 : end do
434 : end do
435 36 : if(tident==1)then
436 36 : identi=isym
437 36 : call wrtout(std_out,sjoin(' symkchk: found identity with number:', itoa(identi)))
438 36 : exit
439 : end if
440 : end do
441 36 : if(tident==0)then
442 0 : errmsg = 'Did not found the identity operation.'
443 0 : ierr = 1; return
444 : end if
445 : end if
446 :
447 : !Here begins the serious business
448 : !The length sorting, etc. of symkpt have been dropped because the
449 : !computational cost is estimated to be negligible.
450 :
451 37 : if(nsym>1 .or. timrev==1)then
452 :
453 : ! Outer loop over kpts
454 5672 : do ikpt=1,nkpt-1
455 :
456 : ! Loop on the symmetries
457 : ! For each k-point and each symmetry transformation, a matching
458 : ! k-point must be found, modulo time reversal if appropriate
459 231509 : do isym=1,nsym
460 :
461 : ! Get the symmetric of the vector
462 903348 : do ii=1,3
463 : ksym(ii)= kptns(1,ikpt)*symrec(ii,1,isym)&
464 : & +kptns(2,ikpt)*symrec(ii,2,isym)&
465 903348 : & +kptns(3,ikpt)*symrec(ii,3,isym)
466 : end do
467 :
468 : ! Second loop k-points
469 133743336 : do ikpt2=1,nkpt
470 :
471 : ! Test for match of symmetric and any vector (including original)
472 : imatch=1
473 534973344 : do ii=1,3
474 401230008 : difk= ksym(ii)-kptns(ii,ikpt2)
475 401230008 : reduce=difk-anint(difk)
476 534973344 : if(abs(reduce)>tol8)imatch=0
477 : end do
478 133743336 : if(imatch==1)exit
479 :
480 : ! Test for match with time reversal
481 133855570 : if(timrev==1)then
482 : imatch=1
483 534518932 : do ii=1,3
484 400889199 : difk= ksym(ii)+kptns(ii,ikpt2)
485 400889199 : reduce=difk-anint(difk)
486 534518932 : if(abs(reduce)>tol8)imatch=0
487 : end do
488 133629733 : if(imatch==1)exit
489 : end if
490 :
491 : end do ! End secondary loop over k-points
492 231472 : if (imatch/=1) then
493 : write(errmsg, '(a,a,a,i0,a,i0,a,a,a,a)' )&
494 0 : 'k-point set must have full space-group symmetry',ch10,&
495 0 : 'there is no match for kpt: ',ikpt,' transformed by symmetry: ',isym,ch10,&
496 0 : 'Action: change kptopt to 2 or 3 and/or change or use shiftk',ch10,&
497 0 : 'shiftk = 0 0 0 is always a safe choice.'
498 0 : ierr = 2; return
499 : end if
500 :
501 : end do ! End loop on isym
502 : end do ! End primary loop over k-points
503 :
504 37 : write(msg,'(a)')' symkchk : k-point set has full space-group symmetry.'
505 111 : call wrtout([std_out, ab_out], msg)
506 : end if
507 :
508 37 : end function symkchk
509 : !!***
510 :
511 : !!****f* m_kpts/kpts_sort
512 : !! NAME
513 : !! kpts_sort
514 : !!
515 : !! FUNCTION
516 : !! Order list of k-points according to their norm.
517 : !!
518 : !! SOURCE
519 :
520 0 : subroutine kpts_sort(gprimd, nkpt, kpts)
521 :
522 : !Arguments ------------------------------------
523 : !scalars
524 : integer,intent(in) :: nkpt
525 : !arrays
526 : real(dp),intent(in) :: gprimd(3, 3)
527 : real(dp),intent(inout) :: kpts(3, nkpt)
528 :
529 : !Local variables-------------------------------
530 : !scalars
531 : integer :: ikpt
532 : !arrays
533 0 : integer,allocatable :: iperm(:)
534 0 : real(dp),allocatable :: knorm2(:), kpts_ord(:,:)
535 : ! *************************************************************************
536 :
537 0 : ABI_MALLOC(knorm2, (nkpt))
538 0 : do ikpt=1,nkpt
539 0 : knorm2(ikpt) = dot_product(kpts(:,ikpt), matmul(gprimd, kpts(:, ikpt)))
540 : end do
541 :
542 0 : ABI_MALLOC(iperm, (nkpt))
543 0 : iperm = [(ikpt, ikpt=1, nkpt)]
544 0 : call sort_dp(nkpt, knorm2, iperm, tol12)
545 0 : ABI_FREE(knorm2)
546 :
547 0 : ABI_MALLOC(kpts_ord, (3, nkpt))
548 0 : do ikpt=1,nkpt
549 0 : kpts_ord(:, ikpt) = kpts(:, iperm(ikpt))
550 : end do
551 0 : kpts = kpts_ord
552 :
553 0 : ABI_FREE(iperm)
554 0 : ABI_FREE(kpts_ord)
555 :
556 0 : end subroutine kpts_sort
557 : !!***
558 :
559 : !!****f* m_kpts/kpts_pack_in_stars
560 : !! NAME
561 : !! kpts_pack_in_stars
562 : !!
563 : !! FUNCTION
564 : !! Pack k-points in stars using kmap symmetry table.
565 : !! Rearrange kpts as well as kmap.
566 : !!
567 : !! SOURCE
568 :
569 8 : subroutine kpts_pack_in_stars(nkpt, kpts, kmap)
570 :
571 : !Arguments ------------------------------------
572 : !scalars
573 : integer,intent(in) :: nkpt
574 : !arrays
575 : real(dp),intent(inout) :: kpts(3, nkpt)
576 : integer,intent(inout) :: kmap(6, nkpt)
577 :
578 : !Local variables-------------------------------
579 : !scalars
580 : integer :: ikpt, seen_ibz, ik_start, ik0, nkibz
581 : integer :: ik_ibz, isym_k, trev_k, tsign, g0_k(3)
582 : logical :: isirr_k
583 : !arrays
584 8 : integer,allocatable :: iperm(:), ibz_ids(:), kmap_ord(:,:), star_pos(:,:)
585 : real(dp) :: swap_kpt(3), swap_kmap(6)
586 8 : real(dp),allocatable :: kpts_ord(:,:)
587 : ! *************************************************************************
588 :
589 : ! Order according to ik_ibz index
590 24 : ABI_MALLOC(ibz_ids, (nkpt))
591 16 : ABI_MALLOC(iperm, (nkpt))
592 2573 : ibz_ids = kmap(1, :)
593 7695 : iperm = [(ikpt, ikpt=1, nkpt)]
594 :
595 8 : call sort_int(nkpt, ibz_ids, iperm)
596 8 : ABI_FREE(ibz_ids)
597 :
598 : ! Rearrange items in _ord arrays.
599 24 : ABI_MALLOC(kpts_ord, (3, nkpt))
600 24 : ABI_MALLOC(kmap_ord, (6, nkpt))
601 2565 : do ikpt=1,nkpt
602 10228 : kpts_ord(:, ikpt) = kpts(:, iperm(ikpt))
603 17907 : kmap_ord(:, ikpt) = kmap(:, iperm(ikpt))
604 : end do
605 :
606 : ! We want each star group to start with the point in the IBZ so an extra shuffle is needed.
607 : ! star_pos stores the beginning of the star group and the position of the base k0 for each star.
608 2565 : nkibz = maxval(kmap(1,:))
609 423 : ABI_ICALLOC(star_pos, (2, nkibz))
610 :
611 : seen_ibz = -1
612 2565 : do ikpt=1,nkpt
613 2557 : ik_ibz = kmap_ord(1, ikpt); isym_k = kmap_ord(2, ikpt)
614 10228 : trev_k = kmap_ord(6, ikpt); g0_k = kmap_ord(3:5, ikpt)
615 8601 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
616 2557 : tsign = 1; if (trev_k == 1) tsign = -1
617 2557 : if (ik_ibz /= seen_ibz) then
618 133 : star_pos(1, ik_ibz) = ikpt
619 133 : seen_ibz = ik_ibz
620 : end if
621 2565 : if (isirr_k) star_pos(2, ik_ibz) = ikpt
622 : end do
623 :
624 : ! Now put ik0 in the ik_start slot if needed.
625 141 : do ik_ibz=1, nkibz
626 133 : ik_start = star_pos(1, ik_ibz)
627 133 : ik0 = star_pos(2, ik_ibz)
628 133 : if (ik_start == ik0) cycle
629 :
630 408 : swap_kpt = kpts_ord(:, ik_start)
631 714 : swap_kmap = kmap_ord(:, ik_start)
632 :
633 408 : kpts_ord(:, ik_start) = kpts_ord(:, ik0)
634 714 : kmap_ord(:, ik_start) = kmap_ord(:, ik0)
635 408 : kpts_ord(:, ik0) = swap_kpt
636 753 : kmap_ord(:, ik0) = swap_kmap
637 : end do
638 :
639 10236 : kpts = kpts_ord
640 17907 : kmap = kmap_ord
641 :
642 8 : ABI_FREE(star_pos)
643 8 : ABI_FREE(iperm)
644 8 : ABI_FREE(kpts_ord)
645 8 : ABI_FREE(kmap_ord)
646 :
647 8 : end subroutine kpts_pack_in_stars
648 : !!***
649 :
650 : !!****f* m_kpts/kpts_map
651 : !! NAME
652 : !! kpts_map
653 : !!
654 : !! FUNCTION
655 : !! Compute symmetry table.
656 : !!
657 : !! INPUTS
658 : !!
659 : !! OUTPUT
660 : !!
661 : !! SOURCE
662 :
663 87045 : integer function kpts_map(mode, kptopt, cryst, krank, nkpt2, kpt2, map, &
664 : qpt, dksqmax_tol) result(ierr) ! optional
665 :
666 : !Arguments ------------------------------------
667 : !scalars
668 : character(len=*),intent(in) :: mode
669 : integer,intent(in) :: kptopt, nkpt2
670 : class(crystal_t),intent(in) :: cryst
671 : class(krank_t),intent(inout) :: krank
672 : real(dp),optional,intent(in) :: dksqmax_tol
673 : !arrays
674 : real(dp),intent(in) :: kpt2(3, nkpt2)
675 : real(dp),optional,intent(in) :: qpt(3)
676 : integer,intent(out) :: map(6, nkpt2)
677 :
678 : !Local variables-------------------------------
679 : !scalars
680 : real(dp) :: dksqmax, my_tol
681 : integer :: timrev, nsym
682 : !arrays
683 : real(dp) :: my_qpt(3)
684 : ! *************************************************************************
685 :
686 87045 : my_qpt = zero; if (present(qpt)) my_qpt = qpt
687 87045 : timrev = kpts_timrev_from_kptopt(kptopt)
688 : ! if no spatial symm. set nsym to 1 to suppress the use of spatial symm.
689 : ! the first symm. is always the identity
690 87045 : if (kptopt==2 .or. kptopt==3) then
691 0 : nsym = 1
692 : else
693 87045 : nsym = cryst%nsym
694 : end if
695 :
696 86232 : select case (mode)
697 : case ("symrel")
698 : ! Note symrel and use_symrec = .False.
699 : ! These are the conventions for the symmetrization of the wavefunctions used in cgtk_rotate.
700 : call krank%get_mapping(nkpt2, kpt2, dksqmax, cryst%gmet, map, &
701 86232 : nsym, cryst%symafm, cryst%symrel, timrev, use_symrec=.False., qpt=my_qpt)
702 :
703 : case ("symrec")
704 : ! Note symrec and use_symrec = .True.
705 : ! These are the conventions for the symmetrization of the DVDB as well as the conventions
706 : ! used in several BZ routines.
707 : ! In principle one should always use this convention but then the symmetrization of the
708 : ! wavefunctions should be rewritten almost completely.
709 :
710 : call krank%get_mapping(nkpt2, kpt2, dksqmax, cryst%gmet, map, &
711 813 : nsym, cryst%symafm, cryst%symrec, timrev, use_symrec=.True., qpt=my_qpt)
712 :
713 : case default
714 87045 : ABI_ERROR(sjoin("Invalid mode:", mode))
715 : end select
716 :
717 87045 : my_tol = tol12; if (present(dksqmax_tol)) my_tol = dksqmax_tol
718 :
719 87045 : ierr = merge(1, 0, dksqmax > my_tol)
720 : !if (ierr /= 0) call wrtout(std_out, sjoin(" CRITICAL WARNING: dksqmax ", ftoa(dksqmax), " > ", ftoa(my_tol)))
721 :
722 : #if 0
723 : !if (present(pack_in_stars)) then
724 : ! if (pack_in_stars)
725 : ! call kpts_pack_in_stars(nkpt2, kpt2, kmap)
726 : ! end if
727 : !end if
728 :
729 :
730 : kmap%nkpt = nkpt2
731 : kmap%kptopt = kptopt
732 : kmap%ierr = ierr
733 : kmap%dksqmax = dksqmax
734 : kmap%mode = trim(mode)
735 : kmap%qpt = qpt
736 : ABI_MALLOC(kmap%inds, (nkpt2))
737 : do ii=1,nkpt2
738 : kmap(ii)%ibz_idx = indkk_kq(1, ii)
739 : kmap(ii)%isym = indkk_kq(2, ii)
740 : kmap(ii)%trev = indkk_kq(6, ii)
741 : kmap(ii)%g0 = indkk_kq(3:5, ii)
742 : kmap(ii)%is_irred = (kmap(ii)%isym == 1 .and. kmap(ii)%trev == 0 .and. all(kmap%g0 == 0))
743 : end do
744 :
745 : !ABI_FREE(kmap%inds)
746 : #endif
747 :
748 87045 : end function kpts_map
749 : !!***
750 :
751 : !!****f* m_kpts/kpts_map_print
752 : !! NAME
753 : !! kpts_map_print
754 : !!
755 : !! FUNCTION
756 : !! Print the symmetry table bz2ibz associated to the BZ bz and the IBZ ibz
757 : !! to a list of units with header. Mode corresponds to the value passed to kpts_map
758 : !! If prtvol is 0, max 20 entries are printed. Use prtvol > 0 to print all k-points.
759 : !!
760 : !! SOURCE
761 :
762 0 : subroutine kpts_map_print(units, header, mode, bz, ibz, bz2ibz, prtvol)
763 :
764 : !Arguments ------------------------------------
765 : !scalars
766 : character(len=*),intent(in) :: header, mode
767 : integer,intent(in) :: prtvol, units(:), bz2ibz(:,:)
768 : real(dp),intent(in) :: bz(:,:), ibz(:,:)
769 :
770 : !Local variables-------------------------------
771 : !scalars
772 : integer :: ik_ibz, ik_bz, isym_k, trev_k, g0_k(3)
773 : logical :: isirr_k
774 : character(len=5000) :: msg
775 : ! *************************************************************************
776 :
777 0 : call wrtout(units, " "//trim(header))
778 0 : select case (mode)
779 : case ("symrec")
780 : call wrtout(units, &
781 0 : " Legend: bz = TS(ibz) + g0 where isym is the index of the symrec operation S and itim is 1 if TR is used.")
782 : case ("symrel")
783 : call wrtout(units, &
784 0 : " Legend: bz = TS^t(ibz) + g0 where isym is the index of the symrel operation S and itim is 1 if TR is used.")
785 : case default
786 0 : ABI_ERROR(sjoin("Invalid mode:", mode))
787 : end select
788 :
789 : ! yes, I'm a barbarian but Fortran string formatting is a pain.
790 0 : msg = " BZ IBZ ibz isym itim g0"
791 0 : call wrtout(units, msg)
792 :
793 0 : do ik_bz=1,size(bz2ibz, dim=2)
794 0 : if (prtvol == 0 .and. ik_bz > 20) then
795 0 : call wrtout(units, "prtvol = 0, max 20 points are written"); exit
796 : end if
797 0 : ik_ibz = bz2ibz(1, ik_bz); isym_k = bz2ibz(2, ik_bz)
798 0 : trev_k = bz2ibz(6, ik_bz); g0_k = bz2ibz(3:5, ik_bz)
799 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
800 : write(msg, '(i6, 2x, 2(a,2x), 3(i4,2x), a)' ) &
801 0 : ik_bz, trim(ktoa(bz(:, ik_bz))), trim(ktoa(ibz(:,ik_ibz))), ik_ibz, isym_k, trev_k, trim(ltoa(g0_k))
802 0 : call wrtout(units, msg)
803 : end do
804 :
805 0 : call wrtout(units, ch10)
806 :
807 0 : end subroutine kpts_map_print
808 : !!***
809 :
810 : !!****f* m_kpts/listkk
811 : !! NAME
812 : !! listkk
813 : !!
814 : !! FUNCTION
815 : !! Given a list of nkpt1 initial k points kptns1 and a list of nkpt2
816 : !! final k points kptns2, associates each final kpt with a "closest"
817 : !! initial k point (or symmetric thereof, also taking possible umklapp)
818 : !! as determined by a metric gmet, that commutes with the symmetry operations.
819 : !! The algorithm does not scale as nkpt1 times nkpt2, thanks
820 : !! to the ordering of the kptns1 and kptns2 vectors according to their
821 : !! lengths, and comparison first between vectors of similar lengths.
822 : !! Returns indirect indexing list indkk.
823 : !!
824 : !! INPUTS
825 : !! gmet(3,3)=reciprocal space metric (bohr^-2)
826 : !! kptns1(3,nkpt1)=list of initial k points (reduced coordinates)
827 : !! kptns2(3,nkpt2)=list of final k points
828 : !! nkpt1=number of initial k points
829 : !! nkpt2=number of final k points
830 : !! nsym=number of symmetry elements in space group
831 : !! sppoldbl=if 1, no spin-polarisation doubling
832 : !! if 2, spin-polarisation doubling using symafm
833 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
834 : !! symmat(3,3,nsym)=symmetry operations (symrel or symrec, depending on value of use_symrec
835 : !! timrev=1 if the use of time-reversal is allowed; 0 otherwise
836 : !! comm=MPI communicator.
837 : !! [use_symrec]: if present and true, symmat assumed to be symrec, otherwise assumed to be symrel (default)
838 : !!
839 : !! OUTPUT
840 : !! dksqmax=maximal value of the norm**2 of the difference between
841 : !! a kpt2 vector and the closest k-point found from the kptns1 set, using symmetries.
842 : !! indkk(nkpt2*sppoldbl,6)=describe k point number of kpt1 that allows to
843 : !! generate wavefunctions closest to given kpt2
844 : !! if sppoldbl=2, use symafm to generate spin down wfs from spin up wfs
845 : !!
846 : !! indkk(:,1)=k point number of kptns1
847 : !! indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
848 : !! (if 0, means no symmetry operation, equivalent to identity )
849 : !! indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
850 : !! to give kpt1b, that is the closest to kpt2.
851 : !! indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
852 : !!
853 : !! NOTES
854 : !! The tolerances tol12 and tol8 aims at giving a machine-independent ordering.
855 : !! (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
856 : !! The tolerance tol12 is used for each component of the k vectors,
857 : !! and for the length of the vectors while the tolerance tol8 is used for
858 : !! the comparison of the squared lengths of the separate vectors.
859 : !!
860 : !! SOURCE
861 :
862 4740 : subroutine listkk(dksqmax, gmet, indkk, kptns1, kptns2, nkpt1, nkpt2, nsym, sppoldbl, symafm, symmat, timrev, comm, &
863 : use_symrec) ! optional
864 :
865 : !Arguments ------------------------------------
866 : !scalars
867 : integer,intent(in) :: nkpt1,nkpt2,nsym,sppoldbl,timrev,comm
868 : real(dp),intent(out) :: dksqmax
869 : logical,optional,intent(in) :: use_symrec
870 : !arrays
871 : integer,intent(in) :: symafm(nsym),symmat(3,3,nsym)
872 : integer,intent(out) :: indkk(nkpt2*sppoldbl,6)
873 : real(dp),intent(in) :: gmet(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
874 :
875 : !Local variables-------------------------------
876 : !scalars
877 : integer,parameter :: usesym=1, limit=1
878 : integer :: nprocs, my_rank, ierr, isk_start, isk_stop
879 : integer :: l3,ig1,ig2,ig3,ii,ikpg1,ikpt1,ikpt2,ikpt2_done, isk
880 : integer :: ilarger,ismaller,itrial
881 : integer :: isppol,isym,itimrev,jkpt1,jsym,jtime
882 : integer :: nsym_used,timrev_used
883 : real(dp) :: dksq,dksqmn,lk2,llarger,ldiff,lsmaller,ltrial,min_l
884 : !real(dp) :: cpu,wall,gflops
885 : character(len=500) :: msg
886 : !arrays
887 : integer :: dkint(3),jdkint(3),k1int(3),k2int(3)
888 4740 : integer, allocatable :: isort(:), tmp_indkk(:,:)
889 : real(dp) :: tsec(2)
890 : real(dp) :: dk(3),kpg1(3),kpt1a(3),k1(3),k2(3)
891 : !real(dp) :: kasq,ka(3)
892 4740 : real(dp),allocatable :: lkpg1(:),lkpg1_sorted(:)
893 : ! *************************************************************************
894 :
895 4740 : call timab(1091, 1, tsec)
896 : !call cwtime(cpu, wall, gflops, "start")
897 :
898 4740 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
899 :
900 4740 : if (sppoldbl<1 .or. sppoldbl>2) then
901 0 : write(msg, '(a,i0,a)' )'The value of sppoldbl is: ',sppoldbl,', but it should be either 1 or 2.'
902 0 : ABI_BUG(msg)
903 : end if
904 :
905 : ! When usesym=0, the old way of converting the wavefunctions (without using the symmetries), is recovered.
906 4740 : nsym_used=nsym
907 4740 : timrev_used=timrev
908 : if(usesym==0)nsym_used=1
909 : if(usesym==0)timrev_used=0
910 :
911 : ! Precompute the length of the kpt1 vectors, also taking into account possible umklapp vectors
912 4740 : l3 = (2*limit+1)**3
913 2902842 : ABI_CALLOC(lkpg1, (l3*nkpt1))
914 2898102 : ABI_CALLOC(lkpg1_sorted, (l3*nkpt1))
915 14220 : ABI_MALLOC(isort, (l3*nkpt1))
916 2893362 : isort = 0
917 :
918 4740 : call xmpi_split_work(nkpt1, comm, isk_start, isk_stop)
919 : !write(std_out,*)' List of kpt1 vectors'; write(std_out,*)' Length of the kpt1 vectors:'
920 :
921 : !$OMP PARALLEL DO PRIVATE(k1, k1int, kpg1, ikpg1)
922 110078 : do ikpt1=isk_start,isk_stop
923 421352 : k1(:) = kptns1(:,ikpt1) !; write(std_out,*)ikpt1,k1(:)
924 421352 : k1int(:) = nint(k1(:) + tol12)
925 421352 : k1(:) = k1(:) - k1int(:)
926 426092 : do ig3=-limit,limit
927 316014 : kpg1(3) = k1(3) + ig3
928 1369394 : do ig2=-limit,limit
929 948042 : kpg1(2) = k1(2) + ig2
930 4108182 : do ig1=-limit,limit
931 2844126 : kpg1(1) = k1(1) + ig1
932 :
933 2844126 : ikpg1 = ig1 + limit + 1 + (2*limit+1)*(ig2+limit) + (2*limit+1)**2*(ig3+limit) + l3*(ikpt1-1)
934 : ! Compute the norm of the vector (also taking into account possible umklapp)
935 : lkpg1(ikpg1) = sqrt(gmet(1,1)*kpg1(1)**2+gmet(2,2)*kpg1(2)**2 + &
936 : gmet(3,3)*kpg1(3)**2+two*(gmet(2,1)*kpg1(2)*kpg1(1) + &
937 2844126 : gmet(3,2)*kpg1(3)*kpg1(2)+gmet(3,1)*kpg1(3)*kpg1(1)))
938 2844126 : lkpg1_sorted(ikpg1) = lkpg1(ikpg1)
939 3792168 : isort(ikpg1) = ikpg1
940 : !write(std_out,*)' ikpt1,ig1,ig2,ig3,lkpg1=',ikpt1,ig1,ig2,ig3,lkpg1(ikpg1)
941 : end do
942 : end do
943 : end do
944 : end do
945 :
946 4740 : if (nprocs > 1) then
947 84 : call xmpi_sum(lkpg1_sorted, comm, ierr)
948 84 : call xmpi_sum(lkpg1, comm, ierr)
949 84 : call xmpi_sum(isort, comm, ierr)
950 : end if
951 : !call cwtime_report(" listkk_loop1", cpu, wall, gflops)
952 :
953 4740 : call sort_dp(l3*nkpt1, lkpg1_sorted, isort, tol12)
954 : ! From "precompute" to "sort_dp" represents more than 50% of the overall wall time for large meshes.
955 : !call cwtime_report(" listkk_sort", cpu, wall, gflops)
956 :
957 : !write(std_out,*)' listkk : output list of kpt1 for checking purposes '
958 : !write(std_out,*)' ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii)) '
959 : !do ii=1,l3*nkpt1
960 : ! ikpt1=(isort(ii)-1)/l3+1
961 : ! write(std_out,*)ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii))
962 : !enddo
963 :
964 4740 : dksqmax = zero
965 2059560 : indkk = 0
966 : ! TODO: Should change API to use this shape.
967 : ! workspace array for improved memory access.
968 14220 : ABI_MALLOC(tmp_indkk, (6, nkpt2*sppoldbl))
969 2368850 : tmp_indkk = 0
970 :
971 : ! Split loop in contiguous blocks
972 4740 : call xmpi_split_work(sppoldbl * nkpt2, comm, isk_start, isk_stop)
973 :
974 9480 : do isppol=1,sppoldbl
975 347210 : do ikpt2=1,nkpt2
976 337730 : isk = ikpt2 + (isppol-1)*nkpt2
977 337730 : if (isk < isk_start .or. isk > isk_stop) cycle
978 :
979 319458 : ikpt2_done=0
980 : ! Precompute the length of the kpt2 vector, with the Umklapp vector such that it is the closest to the Gamma point
981 1277832 : k2(:)=kptns2(:,ikpt2)
982 1277832 : k2int(:)=nint(k2(:)+tol12)
983 1277832 : k2(:)=k2(:)-k2int(:)
984 : lk2=sqrt(gmet(1,1)*k2(1)**2+gmet(2,2)*k2(2)**2+&
985 : gmet(3,3)*k2(3)**2+two*(gmet(2,1)*k2(2)*k2(1)+&
986 319458 : gmet(3,2)*k2(3)*k2(2)+gmet(3,1)*k2(3)*k2(1)))
987 : ! write(std_out, '(a,i4,7es16.6)' )' listkk : ikpt2,kptns2(:,ikpt2),k2(:),lk2=',ikpt2,kptns2(:,ikpt2),k2(:),lk2
988 :
989 : ! Find the kpt1 vector whose length is the most similar to the length of lk2 up to a tolerance.
990 : ! Use a bisection algorithm.
991 319458 : ismaller=0; lsmaller=zero
992 319458 : ilarger=l3*nkpt1+1; llarger=huge(one)
993 :
994 : ! This loop should never reach l3*nkpt1, since this is a bisection algorithm
995 2677339 : do ii=1,l3*nkpt1
996 2677339 : if((ilarger-ismaller)<2 .or. (llarger-lsmaller)<2*tol12)exit
997 2357881 : itrial=(ilarger+ismaller)/2 ; ltrial=lkpg1_sorted(itrial)
998 2677339 : if((ltrial-lk2)>tol12)then
999 : ilarger=itrial ; llarger=ltrial
1000 673016 : else if((ltrial-lk2)<-tol12)then
1001 : ismaller=itrial ; lsmaller=ltrial
1002 : else
1003 316836 : ismaller=itrial ; lsmaller=ltrial
1004 316836 : ilarger=itrial ; llarger=ltrial
1005 : end if
1006 : end do
1007 319458 : itrial=ismaller
1008 319458 : if(abs(llarger-lk2)<abs(lsmaller-lk2)-tol12)itrial=ilarger
1009 319458 : if(itrial==0)itrial=ilarger
1010 319458 : ismaller=itrial ; ilarger=itrial
1011 : !write(std_out,*)' listkk : starting search at itrial=',itrial
1012 :
1013 319458 : dksqmn=huge(one)
1014 :
1015 : ! The ii index is dummy. This avoids an infinite loop.
1016 1775869 : do ii=1,l3*nkpt1
1017 : ! If the difference in length between the trial vector and the target vector is bigger
1018 : ! than the already achieved distance, the search is finished ...
1019 1775868 : ldiff = abs(lkpg1_sorted(itrial) - lk2)
1020 : ! write(std_out,*)' listkk : ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,&
1021 : ! dksqmn=',ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,dksqmn
1022 :
1023 1775868 : if (ldiff**2 > dksqmn+tol8) exit
1024 :
1025 : ! If this k-point has already been examined in a previous batch, skip it
1026 : ! First, compute the minimum of the difference of length of the sets of
1027 : ! associated vectors thanks to Umklapp vectors with the target vector
1028 1550473 : ikpt1 = (isort(itrial)-1) /l3 + 1
1029 44963717 : min_l = minval(abs(lkpg1((ikpt1-1)*l3+1:(ikpt1-1)*l3+l3)-lk2))
1030 :
1031 : ! Then compare with the current ldiff
1032 : ! write(std_out,*)' listkk : ikpt1,min_l,ldiff=',ikpt1,min_l,ldiff
1033 1550473 : if (min_l > ldiff-tol12) then
1034 :
1035 : ! Now, will examine the trial vector, and the symmetric ones
1036 : ! MG FIXME: Here there's a possible problem with the order of symmetries because
1037 : ! in symkpt, time-reversal is the innermost loop. This can create inconsistencies in the symmetry tables.
1038 : ! Besides, one should use symrel^{-1 T} to keep the correspondence between isym -> R or S
1039 4024304 : do itimrev=0,timrev_used
1040 65612902 : do isym=1,nsym_used
1041 :
1042 : ! Select magnetic characteristic of symmetries
1043 63075608 : if (isppol == 1 .and. symafm(isym) == -1) cycle
1044 63057208 : if (isppol == 2 .and. symafm(isym) == 1) cycle
1045 :
1046 : ! Compute symmetric point to kpt1
1047 : if (usesym==1) then
1048 : ! original code only used transpose(symrel)
1049 63057208 : if (present(use_symrec)) then
1050 60085771 : if (use_symrec) then
1051 82845425 : kpt1a(:) = MATMUL(symmat(:,:,isym),kptns1(:,ikpt1))
1052 : else
1053 1419298850 : kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1))
1054 : end if
1055 : else
1056 74285925 : kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1))
1057 : end if
1058 252228832 : kpt1a(:)=(1-2*itimrev)*kpt1a(:)
1059 : else
1060 : kpt1a(:)=kptns1(:,ikpt1)
1061 : end if
1062 :
1063 : ! Compute difference with respect to kpt2, modulo a lattice vector
1064 252228832 : dk(:)=kptns2(:,ikpt2)-kpt1a(:)
1065 : if (usesym==1) then
1066 : ! The tolerance insure similar behaviour on different platforms
1067 : ! XG120418: Actually, *assumes* that the closest point will have reduced
1068 : ! coordinates differing by less than 1/2. There might be elongated cells where this is not correct ...
1069 252228832 : dkint(:)=nint(dk(:)+tol12)
1070 252228832 : dk(:)=dk(:)-dkint(:)
1071 : else
1072 : dkint(:)=0
1073 : end if
1074 :
1075 : ! Compute norm of the difference vector, and update kpt1 if better.
1076 : dksq=gmet(1,1)*dk(1)**2+gmet(2,2)*dk(2)**2+ &
1077 : gmet(3,3)*dk(3)**2+two*(gmet(2,1)*dk(2)*dk(1)+ &
1078 63057208 : gmet(3,2)*dk(3)*dk(2)+gmet(3,1)*dk(3)*dk(1))
1079 :
1080 63057208 : if (dksq < dksqmn+tol8) then
1081 : ! If exactly the right point (without using symmetries neither umklapp vector), will exit the search
1082 : ! Note that in this condition, each coordinate is tested separately, without squaring.
1083 : ! So, it is a much stronger condition than dksqmn < tol12
1084 15208248 : if (sum(abs(kptns2(:,ikpt2)-kptns1(:,ikpt1)))<3*tol12) ikpt2_done = 1
1085 :
1086 : ! Update in three cases: either if succeeded to have exactly the vector, or the distance is better,
1087 : ! or the distance is only slightly worsened so select the lowest itimrev, isym or ikpt1,
1088 : ! in order to respect previous ordering
1089 : if (ikpt2_done==1 .or. &
1090 3708000 : dksq+tol12<dksqmn .or. &
1091 : ( abs(dksq-dksqmn)<tol12 .and. &
1092 : ((itimrev<jtime) .or. &
1093 : (itimrev==jtime .and. isym<jsym) .or. &
1094 : (itimrev==jtime .and. isym==jsym .and. ikpt1<jkpt1))))then
1095 :
1096 1342797 : dksqmn = dksq
1097 1342797 : jkpt1 = ikpt1
1098 1342797 : jsym = isym
1099 1342797 : jtime = itimrev
1100 1342797 : jdkint(:) = dkint(:)
1101 :
1102 : !if (ikpt2_done == 1) then
1103 : ! write(std_out,*)'Succeeded to lower dskmn,ikpt2_done=',dksqmn,ikpt2_done
1104 : ! write(std_out,*)' ikpt1,ikpt2=',ikpt1, ikpt2
1105 : ! write(std_out,*)' ikpt1,isym,dkint(:),itimrev=',ikpt1,isym,dkint(:),itimrev
1106 : ! ka(:) = kpt1a(:) + dkint(:)
1107 : ! kasq=gmet(1,1)*ka(1)**2+gmet(2,2)*ka(2)**2+&
1108 : ! gmet(3,3)*ka(3)**2+two*(gmet(2,1)*ka(2)*ka(1)+&
1109 : ! gmet(3,2)*ka(3)*ka(2)+gmet(3,1)*ka(3)*ka(1))
1110 : ! write(std_out,*)' k1 = ',kpt1a(:)
1111 : ! write(std_out,*)' dkint = ',dkint(:)
1112 : ! write(std_out,*)' Actual k1 = ',ka(:)
1113 : ! write(std_out,*)' k2 = ',kptns2(:,ikpt2)
1114 : ! write(std_out,*)' Actual k1sq = ',kasq
1115 : !end if
1116 : end if
1117 : end if
1118 :
1119 65594502 : if (ikpt2_done==1) exit
1120 : end do ! isym
1121 4024304 : if (ikpt2_done==1) exit
1122 : end do ! itimrev
1123 1487010 : if (ikpt2_done==1) exit
1124 : end if
1125 :
1126 : ! Update the interval that has been explored
1127 1456411 : if (itrial < ismaller) ismaller = itrial
1128 1456411 : if (itrial > ilarger) ilarger = itrial
1129 :
1130 : ! Select the next index to be tried (preferably the smaller indices, but this is a bit arbitrary).
1131 : ! write(std_out,*)' before choosing the next index :'
1132 : ! write(std_out,*)' ismaller,itrial,ilarger=',ismaller,itrial,ilarger
1133 : ! write(std_out,*)' lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1)=',&
1134 : ! lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1)
1135 :
1136 1456411 : if (ismaller>1 .and. ilarger<l3*nkpt1) then
1137 1446557 : if (abs(lkpg1_sorted(ismaller-1)-lk2) < abs(lkpg1_sorted(ilarger+1)-lk2)+tol12) then
1138 : itrial = ismaller-1
1139 : else
1140 760427 : itrial = ilarger+1
1141 : end if
1142 : end if
1143 1456411 : if (ismaller==1 .and. ilarger<l3*nkpt1) itrial = ilarger+1
1144 1775869 : if (ismaller>1 .and. ilarger==l3*nkpt1) itrial = ismaller-1
1145 : !if(ismaller==1 .and. ilarger==l3*nkpt1), we are done with the loop !
1146 : end do ! ikpt1
1147 :
1148 : ! Store indices (lots of cache miss here)
1149 : !indkk(isk, 1) = jkpt1
1150 : !indkk(isk, 2) = jsym
1151 : !indkk(isk, 3:5) = jdkint(:)
1152 : !indkk(isk, 6) = jtime
1153 :
1154 319458 : tmp_indkk(1, isk) = jkpt1
1155 319458 : tmp_indkk(2, isk) = jsym
1156 1277832 : tmp_indkk(3:5, isk) = jdkint(:)
1157 319458 : tmp_indkk(6, isk) = jtime
1158 :
1159 319458 : dksqmax = max(dksqmax, dksqmn)
1160 :
1161 324198 : if (dksqmn < -tol12) then
1162 0 : write(msg, '(a,es16.6)' )'The minimum square of dk has negative norm: dksqmn= ',dksqmn
1163 0 : ABI_BUG(msg)
1164 : end if
1165 :
1166 : ! DEBUG SECTION
1167 : !if (dksqmn > tol5) then
1168 : ! if (present(use_symrec)) then
1169 : ! if (use_symrec) then
1170 : ! kpt1a(:) = MATMUL(symmat(:,:,jsym),kptns1(:,jkpt1))
1171 : ! else
1172 : ! kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,jsym)),kptns1(:,jkpt1))
1173 : ! end if
1174 : ! else
1175 : ! kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,jsym)),kptns1(:,jkpt1))
1176 : ! end if
1177 : ! kpt1a(:)=(1-2*jtime)*kpt1a(:)
1178 : ! print *, "Cannot find k2: ", k2(:)
1179 : ! print *, "Rotated TS(k1): ", kpt1a(:)
1180 : ! print *, "with k1: ", kptns1(:, jkpt1)
1181 : ! print *, "dksqmn: ", dksqmn
1182 : !end if
1183 : !END DEBUG
1184 :
1185 : !write(std_out,'(a,i6,i2,2x,i6,5i3,es24.14)' )' listkk: ikpt2,isppol,indkk(isk,:)=',ikpt2,isppol,indkk(isk,:),dksqmn
1186 : end do ! ikpt2
1187 : end do ! isppol
1188 :
1189 4740 : ABI_FREE(isort)
1190 4740 : ABI_FREE(lkpg1)
1191 4740 : ABI_FREE(lkpg1_sorted)
1192 :
1193 2059560 : indkk = transpose(tmp_indkk)
1194 4740 : ABI_FREE(tmp_indkk)
1195 4740 : if (nprocs > 1) then
1196 84 : call xmpi_sum(indkk, comm, ierr)
1197 84 : dksqmn = dksqmax
1198 84 : call xmpi_max(dksqmn, dksqmax, comm, ierr)
1199 : end if
1200 :
1201 4740 : call timab(1091, 2, tsec)
1202 : !call cwtime_report(" listkk_end", cpu, wall, gflops)
1203 :
1204 4740 : end subroutine listkk
1205 : !!***
1206 :
1207 : !!****f* m_kpts/getkgrid
1208 : !! NAME
1209 : !! getkgrid
1210 : !!
1211 : !! FUNCTION
1212 : !! Compute the grid of k points in the irreducible Brillouin zone.
1213 : !! Note that nkpt (and nkpthf) can be computed by calling this routine with nkpt=0, provided that kptopt/=0.
1214 : !! If downsampling is present, also compute a downsampled k grid.
1215 : !!
1216 : !! INPUTS
1217 : !! chksymbreak= if 1, will check whether the k point grid is symmetric (for kptopt=1,2 and 4), and stop if not.
1218 : !! iout=unit number for echoed output . 0 if no output is wished.
1219 : !! iscf= ( <= 0 =>non-SCF), >0 => SCF) MG: FIXME I don't understand why we have to pass the value iscf.
1220 : !! kptopt=option for the generation of k points. defines whether spatial symmetries and/or time-reversal can be used)
1221 : !! msym=default maximal number of symmetries
1222 : !! nsym=number of symmetries
1223 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1224 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
1225 : !! symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
1226 : !! vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
1227 : !! [downsampling(3) = input variable that governs the downsampling]
1228 : !!
1229 : !! OUTPUT
1230 : !! kptrlen=length of the smallest real space supercell vector associated with the lattice of k points.
1231 : !! nkpt_computed=number of k-points in the IBZ computed in the present routine
1232 : !! If nkpt/=0 the following are also output:
1233 : !! kpt(3,nkpt)=reduced coordinates of k points.
1234 : !! wtk(nkpt)=weight assigned to each k point.
1235 : !! [fullbz(3,nkpt_fullbz)]=k-points generated in the full Brillouin zone.
1236 : !! In output: allocated array with the list of k-points in the BZ.
1237 : !! [kpthf(3,nkpthf)]=k-points generated in the full Brillouin zone, possibly downsampled (for Fock).
1238 : !!
1239 : !! NOTES
1240 : !! msym not needed since nsym is the last index.
1241 : !!
1242 : !! SIDE EFFECTS
1243 : !! Input/Output
1244 : !! nkpt=number of k points (might be zero, see output description)
1245 : !! kptrlatt(3,3)=k-point lattice specification
1246 : !! nshiftk=actual number of k-point shifts in shiftk
1247 : !! shiftk(3,MAX_NSHIFTK)=shift vectors for k point generation
1248 : !! [nkpthf] = number of k points in the full BZ, for the Fock operator.
1249 : !!
1250 : !! SOURCE
1251 :
1252 19764 : subroutine getkgrid(chksymbreak,iout,iscf,kpt,kptopt,kptrlatt,kptrlen,&
1253 19764 : & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk,&
1254 9882 : & fullbz,nkpthf,kpthf,downsampling) ! optional
1255 :
1256 : !Arguments ------------------------------------
1257 : !scalars
1258 : integer,intent(in) :: chksymbreak,iout,iscf,kptopt,msym,nkpt,nsym
1259 : integer,intent(inout),optional :: nkpthf
1260 : integer,intent(inout) :: nshiftk
1261 : integer,intent(inout) :: nkpt_computed !vz_i
1262 : real(dp),intent(out) :: kptrlen
1263 : !arrays
1264 : integer,intent(in) :: symafm(msym),symrel(3,3,msym),vacuum(3)
1265 : integer,optional,intent(in) :: downsampling(3)
1266 : integer,intent(inout) :: kptrlatt(3,3)
1267 9882 : integer,allocatable :: indkpt(:)
1268 9882 : integer,allocatable :: bz2ibz_smap(:,:)
1269 : real(dp),intent(in) :: rprimd(3,3)
1270 : real(dp),intent(inout) :: shiftk(3,MAX_NSHIFTK)
1271 : real(dp),intent(inout) :: kpt(3,nkpt) !vz_i
1272 : real(dp),intent(inout) :: wtk(nkpt)
1273 : real(dp),optional,allocatable,intent(out) :: fullbz(:,:)
1274 : real(dp),optional,intent(out) :: kpthf(:,:)
1275 :
1276 : !Local variables-------------------------------
1277 9882 : real(dp),allocatable :: kpt_tmp(:,:), wtk_tmp(:)
1278 :
1279 : call getkgrid_low(chksymbreak,iout,iscf,kpt_tmp,kptopt,kptrlatt,kptrlen,&
1280 : msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk_tmp,indkpt,bz2ibz_smap,&
1281 12498 : fullbz,nkpthf,kpthf,downsampling)
1282 :
1283 9882 : if (nkpt > 0) then
1284 472267 : kpt(:,1:nkpt) = kpt_tmp(:,1:nkpt)
1285 121135 : wtk(1:nkpt) = wtk_tmp(1:nkpt)
1286 : end if
1287 :
1288 9882 : ABI_SFREE(kpt_tmp)
1289 9882 : ABI_SFREE(wtk_tmp)
1290 9882 : ABI_SFREE(indkpt)
1291 9882 : ABI_SFREE(bz2ibz_smap)
1292 :
1293 9882 : end subroutine getkgrid
1294 : !!***
1295 :
1296 : !!****f* m_kpts/getkgrid_low
1297 : !! NAME
1298 : !! getkgrid_low
1299 : !!
1300 : !! FUNCTION
1301 : !! Compute the grid of k points in the irreducible Brillouin zone.
1302 : !! Note that nkpt (and nkpthf) can be computed by calling this routine with nkpt=0, provided that kptopt/=0.
1303 : !! If downsampling is present, also compute a downsampled k grid.
1304 : !!
1305 : !! INPUTS
1306 : !! chksymbreak= if 1, will check whether the k point grid is symmetric (for kptopt=1,2 and 4), and stop if not.
1307 : !! iout=unit number for echoed output . 0 if no output is wished.
1308 : !! iscf= ( <= 0 =>non-SCF), >0 => SCF) MG: FIXME I don't understand why we have to pass the value iscf.
1309 : !! kptopt=option for the generation of k points (defines whether spatial symmetries and/or time-reversal can be used)
1310 : !! msym=default maximal number of symmetries
1311 : !! nsym=number of symmetries
1312 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1313 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
1314 : !! symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
1315 : !! vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
1316 : !! [downsampling(3) = input variable that governs the downsampling]
1317 : !!
1318 : !! OUTPUT
1319 : !! kptrlen=length of the smallest real space supercell vector associated with the lattice of k points.
1320 : !! nkpt_computed=number of k-points in the IBZ computed in the present routine
1321 : !! If nkpt/=0 the following are also output:
1322 : !! kpt(3,nkpt)=reduced coordinates of k points.
1323 : !! wtk(nkpt)=weight assigned to each k point.
1324 : !! bz2ibz_smap(nkbz, 6)= Mapping BZ --> IBZ.
1325 : !! [fullbz(3,nkpt_fullbz)]=k-points generated in the full Brillouin zone.
1326 : !! In output: allocated array with the list of k-points in the BZ.
1327 : !! [kpthf(3,nkpthf)]=k-points generated in the full Brillouin zone, possibly downsampled (for Fock).
1328 : !!
1329 : !! NOTES
1330 : !! msym not needed since nsym is the last index.
1331 : !!
1332 : !! SIDE EFFECTS
1333 : !! Input/Output
1334 : !! nkpt=number of k points (might be zero, see output description)
1335 : !! kptrlatt(3,3)=k-point lattice specification
1336 : !! nshiftk=actual number of k-point shifts in shiftk
1337 : !! shiftk(3,MAX_NSHIFTK)=shift vectors for k point generation
1338 : !! [nkpthf] = number of k points in the full BZ, for the Fock operator.
1339 : !!
1340 : !! SOURCE
1341 :
1342 10411 : subroutine getkgrid_low(chksymbreak,iout,iscf,kpt,kptopt,kptrlatt,kptrlen,&
1343 10411 : & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk,indkpt,bz2ibz_smap,&
1344 10411 : & fullbz,nkpthf,kpthf,downsampling) ! optional
1345 :
1346 : !Arguments ------------------------------------
1347 : !scalars
1348 : integer,intent(in) :: chksymbreak,iout,iscf,kptopt,msym,nkpt,nsym
1349 : integer,intent(inout),optional :: nkpthf
1350 : integer,intent(inout) :: nshiftk
1351 : integer,intent(inout) :: nkpt_computed !vz_i
1352 : real(dp),intent(out) :: kptrlen
1353 : !arrays
1354 : integer,intent(in) :: symafm(msym),symrel(3,3,msym),vacuum(3)
1355 : integer,optional,intent(in) :: downsampling(3)
1356 : integer,intent(inout) :: kptrlatt(3,3)
1357 : real(dp),intent(in) :: rprimd(3,3)
1358 : real(dp),intent(inout) :: shiftk(3,MAX_NSHIFTK)
1359 : integer,allocatable,intent(out) :: indkpt(:)
1360 : integer,allocatable,intent(out) :: bz2ibz_smap(:,:)
1361 : real(dp),allocatable,intent(out) :: kpt(:,:) !vz_i
1362 : real(dp),allocatable,intent(out) :: wtk(:)
1363 : real(dp),optional,allocatable,intent(out) :: fullbz(:,:)
1364 : real(dp),optional,intent(out) :: kpthf(:,:)
1365 :
1366 : !Local variables-------------------------------
1367 : !scalars
1368 : integer, parameter :: max_number_of_prime=47
1369 : integer :: brav,decreased,found,ii,ikpt,iprime,ishiftk,isym,jshiftk,kshiftk,mkpt,mult
1370 : integer :: nkpthf_computed,nkpt_fullbz,nkptlatt,nshiftk2,nsym_used,option
1371 : integer :: test_prime,timrev
1372 : integer :: nkpt_use
1373 : real(dp) :: length2,ucvol,ucvol_super
1374 : character(len=500) :: msg
1375 : !arrays
1376 : integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,&
1377 : & 29,31,37,41,43, 47,53,59,61,67,&
1378 : & 71,73,79,83,89, 97,101,103,107,109,&
1379 : & 113,127,131,137,139, 149,151,157,163,167,&
1380 : & 173,179,181,191,193, 197,199/)
1381 : integer :: kptrlatt2(3,3)
1382 10411 : integer,allocatable :: belong_chain(:),generator(:),number_in_chain(:)
1383 10411 : integer,allocatable :: repetition_factor(:),symrec(:,:,:)
1384 : ! real(dp) :: cart(3,3)
1385 : real(dp) :: dijk(3),delta_dmult(3),dmult(3),fact_vacuum(3),gmet(3,3)
1386 : real(dp) :: gmet_super(3,3),gprimd(3,3),gprimd_super(3,3),klatt2(3,3)
1387 : real(dp) :: klatt3(3,3),kptrlattr(3,3),ktransf(3,3),ktransf_invt(3,3)
1388 : real(dp) :: metmin(3,3),minim(3,3),rmet(3,3),rmet_super(3,3),rprimd_super(3,3)
1389 10411 : real(dp),allocatable :: deltak(:,:),kpt_fullbz(:,:),shiftk2(:,:),shiftk3(:,:),spkpt(:,:),wtk_folded(:),wtk_fullbz(:)
1390 : ! *************************************************************************
1391 :
1392 10411 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1393 :
1394 : !call cwtime(cpu, wall, gflops, "start")
1395 10411 : if (kptopt==1.or.kptopt==4) then
1396 : ! Cannot use antiferromagnetic symmetry operations to decrease the number of k points
1397 : !XG20191123: now, antiferromagnetic symmetry operations can be used to decrease the number of k points for kptopt==4
1398 7644 : nsym_used=0
1399 246731 : do isym=1,nsym
1400 246731 : if(symafm(isym)==1 .or. kptopt==4)nsym_used=nsym_used+1
1401 : end do
1402 22932 : ABI_MALLOC(symrec,(3,3,nsym_used))
1403 7644 : nsym_used=0
1404 246731 : do isym=1,nsym ! Get the symmetry matrices in terms of reciprocal basis
1405 246731 : if(symafm(isym)==1 .or. kptopt==4)then
1406 237515 : nsym_used=nsym_used+1
1407 237515 : call mati3inv(symrel(:,:,isym),symrec(:,:,nsym_used))
1408 : end if
1409 : end do
1410 2767 : else if (kptopt==2) then
1411 : !Use only the time-reversal
1412 756 : nsym_used=1
1413 756 : ABI_MALLOC(symrec,(3,3,1))
1414 9828 : symrec(1:3,1:3,1)=0
1415 3024 : do ii=1,3
1416 3024 : symrec(ii,ii,1)=1
1417 : end do
1418 : end if
1419 :
1420 10411 : kptrlatt2(:,:)=kptrlatt(:,:)
1421 10411 : nshiftk2=nshiftk
1422 10411 : ABI_MALLOC(shiftk2,(3,MAX_NSHIFTK))
1423 10411 : ABI_MALLOC(shiftk3,(3,MAX_NSHIFTK))
1424 8755651 : shiftk2(:,:)=shiftk(:,:)
1425 :
1426 : !Find a primitive k point lattice, if possible, by decreasing the number of shifts.
1427 10411 : if(nshiftk2/=1)then
1428 :
1429 : do
1430 : ! Loop to be repeated if there has been a successful reduction of nshiftk2
1431 16944 : ABI_MALLOC(deltak,(3,nshiftk2))
1432 16944 : ABI_MALLOC(repetition_factor,(nshiftk2))
1433 11296 : ABI_MALLOC(generator,(nshiftk2))
1434 11296 : ABI_MALLOC(belong_chain,(nshiftk2))
1435 11296 : ABI_MALLOC(number_in_chain,(nshiftk2))
1436 :
1437 5648 : decreased=0
1438 22534 : deltak(1,1:nshiftk2)=shiftk2(1,1:nshiftk2)-shiftk2(1,1)
1439 22534 : deltak(2,1:nshiftk2)=shiftk2(2,1:nshiftk2)-shiftk2(2,1)
1440 22534 : deltak(3,1:nshiftk2)=shiftk2(3,1:nshiftk2)-shiftk2(3,1)
1441 73192 : deltak(:,:)=deltak(:,:)-floor(deltak(:,:)+tol8)
1442 :
1443 : ! Identify for each shift, the smallest repetition prime factor that yields a reciprocal lattice vector.
1444 22534 : repetition_factor(:)=0
1445 5648 : repetition_factor(1)=1
1446 16886 : do ishiftk=2,nshiftk2
1447 19722 : do iprime=1,max_number_of_prime
1448 14014 : test_prime=prime_factor(iprime)
1449 56056 : dmult(:)=test_prime*deltak(:,ishiftk)
1450 56116 : if(sum(abs( dmult(:)-nint(dmult(:)) ))<tol8)then
1451 11178 : repetition_factor(ishiftk)=test_prime
1452 11178 : exit
1453 : end if
1454 : end do
1455 : end do
1456 :
1457 : ! Initialize the selection of tentative generators
1458 22534 : generator(:)=1
1459 22534 : do ishiftk=1,nshiftk2
1460 22534 : if(repetition_factor(ishiftk)==0 .or. repetition_factor(ishiftk)==1)generator(ishiftk)=0
1461 : end do
1462 :
1463 : ! Try different shifts as generators, by order of increasing repetition factor,
1464 : ! provided they are equal or bigger than 2
1465 5650 : do iprime=1,max_number_of_prime
1466 5706 : do ishiftk=2,nshiftk2
1467 : ! Note that ishiftk=1 is never a generator. It is the reference starting point.
1468 5704 : if(generator(ishiftk)==1 .and. repetition_factor(ishiftk)==prime_factor(iprime))then
1469 : ! Test the generator : is it indeed closed ?
1470 5648 : if(prime_factor(iprime)/=2)then
1471 4 : do mult=2,prime_factor(iprime)-1
1472 8 : dmult(:)=mult*deltak(:,ishiftk)
1473 6 : found=0
1474 6 : do jshiftk=1,nshiftk2
1475 24 : delta_dmult(:)=deltak(:,jshiftk)-dmult(:)
1476 24 : if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then
1477 : found=1
1478 : exit
1479 : end if
1480 : end do
1481 4 : if(found==0)exit
1482 : end do
1483 2 : if(found==0)generator(ishiftk)=0
1484 : end if
1485 5648 : if(generator(ishiftk)==0)cycle
1486 : else
1487 : cycle
1488 : end if
1489 : ! Now, test whether all k points can be found in all possible chains
1490 22534 : belong_chain(:)=0
1491 22534 : do jshiftk=1,nshiftk2
1492 : ! Initialize a chain starting from a k point not yet in a chain
1493 22534 : if(belong_chain(jshiftk)==0)then
1494 37032 : number_in_chain(:)=0 ! Not a member of the chain (yet)
1495 8442 : number_in_chain(jshiftk)=1 ! The first point in chain
1496 16886 : do mult=1,prime_factor(iprime)-1
1497 33776 : dmult(:)=mult*deltak(:,ishiftk)
1498 8444 : found=0
1499 8690 : do kshiftk=jshiftk+1,nshiftk2
1500 34760 : delta_dmult(:)=deltak(:,kshiftk)-deltak(:,jshiftk)-dmult(:)
1501 34760 : if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then
1502 8444 : found=1
1503 8444 : number_in_chain(kshiftk)=mult+1
1504 : exit
1505 : end if
1506 : end do
1507 8442 : if(found==0)then
1508 0 : generator(ishiftk)=0
1509 0 : exit
1510 : end if
1511 : end do
1512 8442 : if(generator(ishiftk)==1)then
1513 : ! Store the chain
1514 37032 : do kshiftk=1,nshiftk2
1515 37032 : if(number_in_chain(kshiftk)/=0)belong_chain(kshiftk)=number_in_chain(kshiftk)
1516 : end do
1517 : else
1518 : exit
1519 : end if
1520 : end if
1521 : end do
1522 :
1523 5648 : if(generator(ishiftk)==0)cycle
1524 :
1525 : ! For the generator based on ishiftk, all the k points have been found to belong to one chain.
1526 : ! All the initializing k points in the different chains have belong_chain(:)=1 .
1527 : ! They must be kept, and the others thrown away.
1528 5648 : ktransf(:,:)=0.0_dp
1529 5648 : ktransf(1,1)=1.0_dp
1530 5648 : ktransf(2,2)=1.0_dp
1531 5648 : ktransf(3,3)=1.0_dp
1532 : ! Replace one of the unit vectors by the shift vector deltak(:,ishiftk).
1533 : ! However, must pay attention not to make linear combinations.
1534 : ! Also, choose positive sign for first-non-zero value.
1535 5648 : if(abs(deltak(1,ishiftk)-nint(deltak(1,ishiftk)))>tol8)then
1536 11560 : if(deltak(1,ishiftk)>0)ktransf(:,1)= deltak(:,ishiftk)
1537 2890 : if(deltak(1,ishiftk)<0)ktransf(:,1)=-deltak(:,ishiftk)
1538 2758 : else if(abs(deltak(2,ishiftk)-nint(deltak(2,ishiftk)))>tol8)then
1539 11032 : if(deltak(2,ishiftk)>0)ktransf(:,2)= deltak(:,ishiftk)
1540 2758 : if(deltak(2,ishiftk)<0)ktransf(:,2)=-deltak(:,ishiftk)
1541 0 : else if(abs(deltak(3,ishiftk)-nint(deltak(3,ishiftk)))>tol8)then
1542 0 : if(deltak(3,ishiftk)>0)ktransf(:,3)= deltak(:,ishiftk)
1543 0 : if(deltak(3,ishiftk)<0)ktransf(:,3)=-deltak(:,ishiftk)
1544 : end if
1545 : ! Copy the integers to real(dp)
1546 73424 : kptrlattr(:,:)=kptrlatt2(:,:)
1547 : ! Go to reciprocal space
1548 5648 : call matr3inv(kptrlattr,klatt2)
1549 : ! Make the transformation
1550 22592 : do ii=1,3
1551 73424 : klatt3(:,ii)=ktransf(1,ii)*klatt2(:,1)+ktransf(2,ii)*klatt2(:,2)+ktransf(3,ii)*klatt2(:,3)
1552 : end do
1553 : ! Back to real space
1554 5648 : call matr3inv(klatt3,kptrlattr)
1555 : ! real(dp) to integer
1556 73424 : kptrlatt2(:,:)=nint(kptrlattr(:,:))
1557 : ! Prepare the transformation of the shifts
1558 5648 : call matr3inv(ktransf,ktransf_invt)
1559 5648 : decreased=1
1560 5648 : kshiftk=0
1561 22534 : do jshiftk=1,nshiftk2
1562 22534 : if(belong_chain(jshiftk)==1)then
1563 8442 : kshiftk=kshiftk+1
1564 : ! Place the shift with index jshiftk in place of the one in kshiftk,
1565 : ! also transform the shift from the old to the new coordinate system
1566 : shiftk3(:,kshiftk)=ktransf_invt(1,:)*shiftk2(1,jshiftk)+&
1567 : & ktransf_invt(2,:)*shiftk2(2,jshiftk)+&
1568 33768 : & ktransf_invt(3,:)*shiftk2(3,jshiftk)
1569 : end if
1570 : end do
1571 5648 : nshiftk2=nshiftk2/prime_factor(iprime)
1572 39416 : shiftk2(:,1:nshiftk2)=shiftk3(:,1:nshiftk2)-floor(shiftk3(:,1:nshiftk2)+tol8)
1573 5648 : if(kshiftk/=nshiftk2)then
1574 0 : ABI_BUG('The search for a primitive k point lattice contains a bug.')
1575 : end if
1576 :
1577 : ! If this trial shift was successful, must exit the loop on trial ishiftk,
1578 : ! and reinitialize the global loop
1579 2 : if(decreased==1)exit
1580 : end do ! ishiftk
1581 0 : if(decreased==1)exit
1582 : end do ! iprime
1583 :
1584 5648 : ABI_FREE(belong_chain)
1585 5648 : ABI_FREE(deltak)
1586 5648 : ABI_FREE(number_in_chain)
1587 5648 : ABI_FREE(repetition_factor)
1588 5648 : ABI_FREE(generator)
1589 :
1590 5648 : if(decreased==0 .or. nshiftk2==1)exit
1591 :
1592 : end do ! Infinite loop
1593 :
1594 : end if ! End nshiftk being 1 or larger
1595 :
1596 : !Impose shiftk coordinates to be in [0,1[
1597 20822 : do ishiftk=1,nshiftk2
1598 52055 : do ii=1,3
1599 31233 : if(shiftk2(ii,ishiftk)>one-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)-1.0_dp
1600 41644 : if(shiftk2(ii,ishiftk)<-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)+1.0_dp
1601 : end do
1602 : end do
1603 :
1604 : !Compute the number of k points in the G-space unit cell
1605 : nkptlatt=kptrlatt2(1,1)*kptrlatt2(2,2)*kptrlatt2(3,3) &
1606 : & +kptrlatt2(1,2)*kptrlatt2(2,3)*kptrlatt2(3,1) &
1607 : & +kptrlatt2(1,3)*kptrlatt2(2,1)*kptrlatt2(3,2) &
1608 : & -kptrlatt2(1,2)*kptrlatt2(2,1)*kptrlatt2(3,3) &
1609 : & -kptrlatt2(1,3)*kptrlatt2(2,2)*kptrlatt2(3,1) &
1610 10411 : & -kptrlatt2(1,1)*kptrlatt2(2,3)*kptrlatt2(3,2)
1611 :
1612 : !Check whether the number of k points is positive, otherwise, change the handedness of kptrlatt2
1613 10411 : if(nkptlatt<=0)then
1614 : ! write(std_out,*)' getkgrid : nkptlatt is negative !'
1615 100 : kptrlatt2(:,3)=-kptrlatt2(:,3)
1616 25 : nkptlatt=-nkptlatt
1617 50 : do ishiftk=1,nshiftk2
1618 50 : shiftk2(3,ishiftk)=-shiftk2(3,ishiftk)
1619 : end do
1620 : end if
1621 :
1622 : !Determine the smallest supercell R-vector whose contribution
1623 : !is not taken correctly into account in the k point integration.
1624 : !Increase enormously the size of the cell when vacuum is present.
1625 41644 : fact_vacuum(:)=1
1626 10411 : if(vacuum(1)==1)fact_vacuum(1)=1000.0_dp
1627 10411 : if(vacuum(2)==1)fact_vacuum(2)=1000.0_dp
1628 10411 : if(vacuum(3)==1)fact_vacuum(3)=1000.0_dp
1629 41644 : do ii=1,3
1630 : rprimd_super(:,ii)=fact_vacuum(1)*rprimd(:,1)*kptrlatt2(1,ii)+&
1631 : & fact_vacuum(2)*rprimd(:,2)*kptrlatt2(2,ii)+&
1632 135343 : & fact_vacuum(3)*rprimd(:,3)*kptrlatt2(3,ii)
1633 : end do
1634 :
1635 10411 : call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super)
1636 10411 : call smallprim(metmin,minim,rprimd_super)
1637 10411 : length2=min(metmin(1,1),metmin(2,2),metmin(3,3))
1638 10411 : kptrlen=sqrt(length2)
1639 :
1640 : !write(msg,'(a,es16.6)' )' getkgrid : length of smallest supercell vector (bohr)=',kptrlen
1641 : !call wrtout(std_out,msg)
1642 : ! If the number of shifts has been decreased, determine the set of kptrlatt2 vectors
1643 : ! with minimal length (without using fact_vacuum)
1644 : ! It is worth to determine the minimal set of vectors so that the kptrlatt that is output
1645 : ! does not seem screwy, although correct but surprising.
1646 10411 : if(nshiftk/=nshiftk2)then
1647 11536 : do ii=1,3
1648 37492 : rprimd_super(:,ii)=rprimd(:,1)*kptrlatt2(1,ii)+rprimd(:,2)*kptrlatt2(2,ii)+rprimd(:,3)*kptrlatt2(3,ii)
1649 : end do
1650 2884 : call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super)
1651 : ! Shift vectors in cartesian coordinates (reciprocal space)
1652 5768 : do ishiftk=1,nshiftk2
1653 : shiftk3(:,ishiftk)=gprimd_super(:,1)*shiftk2(1,ishiftk)+&
1654 : & gprimd_super(:,2)*shiftk2(2,ishiftk)+&
1655 14420 : & gprimd_super(:,3)*shiftk2(3,ishiftk)
1656 : end do
1657 2884 : call smallprim(metmin,minim,rprimd_super)
1658 2884 : call metric(gmet_super,gprimd_super,-1,rmet_super,minim,ucvol_super)
1659 : ! This is the new kptrlatt2
1660 11536 : do ii=1,3
1661 : dijk(:)=gprimd(1,:)*minim(1,ii)+&
1662 : & gprimd(2,:)*minim(2,ii)+&
1663 34608 : & gprimd(3,:)*minim(3,ii)
1664 37492 : kptrlatt2(:,ii)=nint(dijk(:))
1665 : end do
1666 : ! Shifts in the new set of kptrlatt vectors
1667 5768 : do ishiftk=1,nshiftk2
1668 : shiftk2(:,ishiftk)=minim(1,:)*shiftk3(1,ishiftk)+&
1669 : & minim(2,:)*shiftk3(2,ishiftk)+&
1670 14420 : & minim(3,:)*shiftk3(3,ishiftk)
1671 : end do
1672 : end if
1673 :
1674 : !brav=1 is able to treat all bravais lattices.
1675 10411 : brav=1
1676 10411 : mkpt=nkptlatt*nshiftk2
1677 :
1678 31233 : ABI_MALLOC(spkpt,(3,mkpt))
1679 10411 : option=0
1680 10411 : if(iout/=0)option=1
1681 :
1682 : !call cwtime_report(' shifts', cpu, wall, gflops)
1683 :
1684 10411 : if (present(downsampling))then
1685 7266 : call smpbz(brav,iout,kptrlatt2,mkpt,nkpthf_computed,nshiftk2,option,shiftk2,spkpt,downsampling=downsampling)
1686 7266 : if (present(kpthf) .and. nkpthf/=0) then
1687 : ! Returns list of k-points in the Full BZ, possibly downsampled for Fock
1688 3006151 : kpthf = spkpt(:,1:nkpthf)
1689 : end if
1690 7266 : nkpthf=nkpthf_computed
1691 : end if
1692 :
1693 10411 : call smpbz(brav,iout,kptrlatt2,mkpt,nkpt_fullbz,nshiftk2,option,shiftk2,spkpt)
1694 : !call cwtime_report(' smpbz', cpu, wall, gflops)
1695 :
1696 10411 : if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then
1697 :
1698 25200 : ABI_MALLOC(indkpt,(nkpt_fullbz))
1699 25200 : ABI_MALLOC(kpt_fullbz,(3,nkpt_fullbz))
1700 25200 : ABI_MALLOC(bz2ibz_smap, (6, nkpt_fullbz))
1701 : #if 1
1702 25200 : ABI_MALLOC(wtk_fullbz,(nkpt_fullbz))
1703 16800 : ABI_MALLOC(wtk_folded,(nkpt_fullbz))
1704 :
1705 21639772 : kpt_fullbz(:,:)=spkpt(:,1:nkpt_fullbz)
1706 5416243 : wtk_fullbz(1:nkpt_fullbz)=1.0_dp/dble(nkpt_fullbz)
1707 :
1708 8400 : timrev=1;if (kptopt==4) timrev=0
1709 :
1710 5416243 : indkpt = 0
1711 : call symkpt(chksymbreak,gmet,indkpt,iout,kpt_fullbz,nkpt_fullbz,&
1712 8400 : & nkpt_computed,nsym_used,symrec,timrev,wtk_fullbz,wtk_folded,bz2ibz_smap,xmpi_comm_self)
1713 :
1714 8400 : ABI_FREE(symrec)
1715 8400 : ABI_FREE(wtk_fullbz)
1716 :
1717 : !do ikpt=1,nkpt_fullbz
1718 : ! write(*,*) ikpt, indkpt(ikpt), bz2ibz_smap(1,ikpt), indkpt(bz2ibz_smap(1,ikpt))
1719 : !end do
1720 : #else
1721 : kpt_fullbz(:,:)=spkpt(:,1:nkpt_fullbz)
1722 :
1723 : timrev=1;if (kptopt==4) timrev=0
1724 :
1725 : call symkpt_new(chksymbreak,gmet,indkpt,iout,kpt_fullbz,nkpt_fullbz,&
1726 : & nkpt_computed,nsym_used,symrec,timrev,bz2ibz_smap,xmpi_comm_self)
1727 :
1728 : ABI_FREE(symrec)
1729 : ABI_CALLOC(wtk_folded,(nkpt_fullbz))
1730 : do ii=1,nkpt_fullbz
1731 : ikpt = indkpt(bz2ibz_smap(1,ii))
1732 : wtk_folded(ikpt) = wtk_folded(ikpt) + one
1733 : end do
1734 : wtk_folded = wtk_folded / nkpt_fullbz
1735 : #endif
1736 :
1737 2011 : else if(kptopt==3)then
1738 686111 : ABI_ICALLOC(bz2ibz_smap, (6, nkpt_fullbz))
1739 196319 : bz2ibz_smap(1,:) = [(ii,ii=1,nkpt_fullbz)]
1740 99165 : bz2ibz_smap(2,:) = 1 !isym
1741 2011 : nkpt_computed=nkpt_fullbz
1742 : end if
1743 : !call cwtime_report(' symkpt', cpu, wall, gflops)
1744 :
1745 : !The number of k points has been computed from kptopt, kptrlatt, nshiftk, shiftk,
1746 : !and the eventual symmetries, it is presently called nkpt_computed.
1747 10411 : nkpt_use = nkpt
1748 10411 : if (nkpt<0) nkpt_use = nkpt_computed
1749 :
1750 : !Check that the argument nkpt is coherent with nkpt_computed, if nkpt/=0.
1751 10411 : if(nkpt_use/=nkpt_computed .and. nkpt/=0)then
1752 : write(msg, '(a,i0,5a,i0,7a)') &
1753 0 : 'The argument nkpt = ',nkpt_use,', does not match',ch10,&
1754 0 : 'the number of k points generated by kptopt, kptrlatt, shiftk,',ch10,&
1755 0 : 'and the eventual symmetries, that is, nkpt= ',nkpt_computed,'.',ch10,&
1756 0 : 'However, note that it might be due to the user,',ch10,&
1757 0 : 'if nkpt is explicitly defined in the input file.',ch10,&
1758 0 : 'In this case, please check your input file.'
1759 0 : ABI_BUG(msg)
1760 : end if
1761 :
1762 31233 : ABI_MALLOC(kpt,(3,nkpt_use))
1763 31233 : ABI_MALLOC(wtk,(nkpt_use))
1764 :
1765 10411 : if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then
1766 :
1767 8400 : if(nkpt_use/=0)then
1768 317274 : do ikpt=1,nkpt_use
1769 1254972 : kpt(:,ikpt)=kpt_fullbz(:,indkpt(ikpt))
1770 317274 : if(iscf>=0 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(ikpt)=wtk_folded(indkpt(ikpt))
1771 : end do
1772 : end if
1773 :
1774 8400 : if (present(fullbz)) then
1775 : ! Returns list of k-points in the Full BZ.
1776 538 : ABI_MOVE_ALLOC(kpt_fullbz,fullbz)
1777 : else
1778 7862 : ABI_FREE(kpt_fullbz)
1779 : end if
1780 :
1781 8400 : ABI_FREE(wtk_folded)
1782 :
1783 2011 : else if(kptopt==3)then
1784 :
1785 2011 : if(nkpt_use/=0)then
1786 224529 : kpt(:,1:nkpt_use)=spkpt(:,1:nkpt_use)
1787 56693 : if(iscf>1 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(1:nkpt_use)=1.0_dp/dble(nkpt_use)
1788 : end if
1789 :
1790 2011 : if (present(fullbz)) then
1791 : ! Returns list of k-points in the Full BZ.
1792 39 : ABI_MALLOC(fullbz,(3,nkpt_fullbz))
1793 3610 : fullbz = spkpt(:,1:nkpt_fullbz)
1794 : end if
1795 :
1796 : end if
1797 :
1798 10411 : ABI_FREE(spkpt)
1799 10411 : kptrlatt(:,:)=kptrlatt2(:,:)
1800 10411 : nshiftk=nshiftk2
1801 52055 : shiftk(:,1:nshiftk)=shiftk2(:,1:nshiftk)
1802 10411 : ABI_FREE(shiftk2)
1803 10411 : ABI_FREE(shiftk3)
1804 :
1805 20822 : end subroutine getkgrid_low
1806 : !!***
1807 :
1808 : !!****f* m_kpts/get_full_kgrid
1809 : !! NAME
1810 : !! get_full_kgrid
1811 : !!
1812 : !! FUNCTION
1813 : !! Create full grid of kpoints and find equivalent
1814 : !! irred ones. Duplicates work in getkgrid, but need all outputs of kpt_fullbz, and indkpt
1815 : !!
1816 : !! INPUTS
1817 : !! kpt(3,nkpt)=irreducible kpoints
1818 : !! kptrlatt(3,3)=lattice vectors for full kpoint grid
1819 : !! nkpt=number of irreducible kpoints
1820 : !! nkpt_fullbz=number of kpoints in full brillouin zone
1821 : !! nshiftk=number of kpoint grid shifts
1822 : !! nsym=number of symmetries
1823 : !! shiftk(3,nshiftk)=kpoint shifts
1824 : !! symrel(3,3,nsym)=symmetry matrices in real space
1825 : !!
1826 : !! OUTPUT
1827 : !! indkpt(nkpt_fullbz)=non-symmetrized indices of the k-points (see symkpt.f)
1828 : !! kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone
1829 : !!
1830 : !! NOTES
1831 : !! MG: The present implementation always assumes kptopt==1 !!!!
1832 : !!
1833 : !! TODO: This routine should be removed
1834 : !!
1835 : !! SOURCE
1836 :
1837 0 : subroutine get_full_kgrid(indkpt,kpt,kpt_fullbz,kptrlatt,nkpt,nkpt_fullbz,nshiftk,nsym,shiftk,symrel)
1838 :
1839 : !Arguments ------------------------------------
1840 : !scalars
1841 : integer,intent(in) :: nkpt,nkpt_fullbz,nshiftk,nsym
1842 : !arrays
1843 : integer,intent(in) :: kptrlatt(3,3),symrel(3,3,nsym)
1844 : integer,intent(out) :: indkpt(nkpt_fullbz)
1845 : real(dp),intent(in) :: kpt(3,nkpt),shiftk(3,nshiftk)
1846 : real(dp),intent(out) :: kpt_fullbz(3,nkpt_fullbz)
1847 :
1848 : !Local variables-------------------------------
1849 : !scalars
1850 : integer :: ikpt,isym,itim,timrev
1851 : integer :: symrankkpt
1852 : character(len=500) :: msg
1853 0 : type(krank_t) :: krank
1854 : !arrays
1855 0 : integer :: inv_symrel(3,3,nsym)
1856 : real(dp) :: k2(3)
1857 : ! *********************************************************************
1858 :
1859 : !Invert symrels => gives symrels for kpoints
1860 :
1861 0 : do isym=1,nsym
1862 0 : call mati3inv (symrel(:,:,isym),inv_symrel(:,:,isym))
1863 : end do
1864 :
1865 0 : call get_kpt_fullbz(kpt_fullbz,kptrlatt,nkpt_fullbz,nshiftk,shiftk)
1866 :
1867 : ! make full k-point rank arrays
1868 0 : call krank%init(nkpt, kpt)
1869 :
1870 : !find equivalence to irred kpoints in kpt
1871 0 : indkpt(:) = 0
1872 0 : timrev=1 ! includes the time inversion symmetry
1873 0 : do ikpt=1,nkpt_fullbz
1874 0 : do isym=1,nsym
1875 0 : do itim=1,(1-2*timrev),-2
1876 :
1877 : k2(:) = itim*(inv_symrel(:,1,isym)*kpt_fullbz(1,ikpt) + &
1878 : inv_symrel(:,2,isym)*kpt_fullbz(2,ikpt) + &
1879 0 : inv_symrel(:,3,isym)*kpt_fullbz(3,ikpt))
1880 :
1881 0 : symrankkpt = krank%get_rank(k2)
1882 0 : if (krank%invrank(symrankkpt) /= -1) indkpt(ikpt) = krank%invrank(symrankkpt)
1883 :
1884 : end do ! loop time reversal symmetry
1885 : end do ! loop sym ops
1886 :
1887 0 : if (indkpt(ikpt) == 0) then
1888 0 : write(msg,'(a,i0)')' indkpt(ikpt) is still 0: no irred kpoint is equiv to ikpt ',ikpt
1889 0 : ABI_BUG(msg)
1890 : end if
1891 : end do ! loop full kpts
1892 :
1893 0 : call krank%free()
1894 :
1895 0 : end subroutine get_full_kgrid
1896 : !!***
1897 :
1898 : !!****f* m_kpts/get_kpt_fullbz
1899 : !! NAME
1900 : !! get_kpt_fullbz
1901 : !!
1902 : !! FUNCTION
1903 : !! Create full grid of kpoints from kptrlatt and shiftk
1904 : !!
1905 : !! INPUTS
1906 : !! kptrlatt(3,3)=lattice vectors for full kpoint grid
1907 : !! nkpt_fullbz=number of kpoints in full brillouin zone
1908 : !! nshiftk=number of kpoint grid shifts
1909 : !! shiftk(3,nshiftk)=kpoint shifts
1910 : !!
1911 : !! OUTPUT
1912 : !! kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone
1913 : !!
1914 : !! SOURCE
1915 :
1916 0 : subroutine get_kpt_fullbz(kpt_fullbz,kptrlatt,nkpt_fullbz,nshiftk,shiftk)
1917 :
1918 : !Arguments ------------------------------------
1919 : !scalars
1920 : integer,intent(in) :: nkpt_fullbz,nshiftk
1921 : !arrays
1922 : integer,intent(in) :: kptrlatt(3,3)
1923 : real(dp),intent(in) :: shiftk(3,nshiftk)
1924 : real(dp),intent(out) :: kpt_fullbz(3,nkpt_fullbz)
1925 :
1926 : !Local variables-------------------------------
1927 : !scalars
1928 : integer, parameter :: max_number_of_prime=47
1929 : integer :: det,ii,ikshft,iprim,jj,kk,nn
1930 : character(len=500) :: msg
1931 : !arrays
1932 : integer :: boundmax(3),boundmin(3),common_factor(3)
1933 : integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,&
1934 : & 29,31,37,41,43, 47,53,59,61,67,&
1935 : & 71,73,79,83,89, 97,101,103,107,109,&
1936 : & 113,127,131,137,139, 149,151,157,163,167,&
1937 : & 173,179,181,191,193, 197,199/)
1938 : real(dp) :: k1(3),k2(3),klatt(3,3),rlatt(3,3),shift(3),test_rlatt(3,3)
1939 : ! *********************************************************************
1940 :
1941 : !Identify first factors that can be used to rescale the three kptrlatt vectors
1942 : !Only test a large set of prime factors, though ...
1943 0 : do jj=1,3
1944 0 : common_factor(jj)=1
1945 0 : rlatt(:,jj)=kptrlatt(:,jj)
1946 0 : do iprim=1,max_number_of_prime
1947 0 : test_rlatt(:,jj)=rlatt(:,jj)/dble(prime_factor(iprim))
1948 : ! If one of the components is lower than 1 in absolute value, then it is not worth to continue the search.
1949 0 : if(minval(abs(abs(test_rlatt(:,jj))-half))<half-tol8)exit
1950 0 : do
1951 0 : if(sum(abs(test_rlatt(:,jj)-nint(test_rlatt(:,jj)) ))<tol8)then
1952 0 : common_factor(jj)=prime_factor(iprim)*common_factor(jj)
1953 0 : rlatt(:,jj)=rlatt(:,jj)/dble(prime_factor(iprim))
1954 0 : test_rlatt(:,jj)=test_rlatt(:,jj)/dble(prime_factor(iprim))
1955 : else
1956 : exit
1957 : end if
1958 : end do
1959 : end do
1960 : end do
1961 0 : call mati3det(kptrlatt,det)
1962 0 : det=det/(common_factor(1)*common_factor(2)*common_factor(3))
1963 :
1964 0 : rlatt(:,:)=kptrlatt(:,:)
1965 0 : call matr3inv(rlatt,klatt)
1966 : !Now, klatt contains the three primitive vectors of the k lattice,
1967 : !in reduced coordinates. One builds all k vectors that
1968 : !are contained in the first Brillouin zone, with coordinates
1969 : !in the interval [0,1[ . First generate boundaries of a big box.
1970 : !In order to generate all possible vectors in the reciprocal space,
1971 : !one must consider all multiples of the primitive ones, until a vector with only integers is found.
1972 : !The maximum bound is the scale of the corresponding kptrlatt vector, times the determinant of kptrlatt. Also consider negative vectors.
1973 : !On this basis, compute the bounds.
1974 0 : do jj=1,3
1975 : ! To accommodate the shifts, boundmin starts from -1
1976 : ! Well, this is not a complete solution ...
1977 0 : boundmin(jj)=-1-common_factor(jj)*abs(det)
1978 0 : boundmax(jj)=common_factor(jj)*abs(det)
1979 : end do
1980 :
1981 0 : nn=1
1982 0 : do kk=boundmin(3),boundmax(3)
1983 0 : do jj=boundmin(2),boundmax(2)
1984 0 : do ii=boundmin(1),boundmax(1)
1985 0 : do ikshft=1,nshiftk
1986 :
1987 : ! Coordinates of the trial k point with respect to the k primitive lattice
1988 0 : k1(1)=ii+shiftk(1,ikshft)
1989 0 : k1(2)=jj+shiftk(2,ikshft)
1990 0 : k1(3)=kk+shiftk(3,ikshft)
1991 :
1992 : ! Reduced coordinates of the trial k point
1993 0 : k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3)
1994 :
1995 : ! Eliminate the point if outside [0,1[
1996 0 : if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle
1997 0 : if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle
1998 0 : if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle
1999 :
2000 : ! Wrap the trial values in the interval ]-1/2,1/2] .
2001 0 : call wrap2_pmhalf(k2(1),k1(1),shift(1))
2002 0 : call wrap2_pmhalf(k2(2),k1(2),shift(2))
2003 0 : call wrap2_pmhalf(k2(3),k1(3),shift(3))
2004 0 : if(nn > nkpt_fullbz) then
2005 0 : write (msg,'(a,i0)')' nkpt_fullbz mis-estimated, exceed nn=',nn
2006 0 : ABI_BUG(msg)
2007 : end if
2008 0 : kpt_fullbz(:,nn)=k1(:)
2009 0 : nn=nn+1
2010 : end do
2011 : end do
2012 : end do
2013 : end do
2014 0 : nn = nn-1
2015 :
2016 0 : if (nn /= nkpt_fullbz) then
2017 0 : write (msg,'(2(a,i0),a,a)')' nkpt_fullbz= ',nkpt_fullbz,' underestimated nn=',nn,&
2018 0 : & ch10, "Perhaps your k grid or shifts do not correspond to the symmetry?"
2019 0 : ABI_BUG(msg)
2020 : end if
2021 :
2022 0 : end subroutine get_kpt_fullbz
2023 : !!***
2024 :
2025 : !!****f* m_kpts/smpbz
2026 : !! NAME
2027 : !! smpbz
2028 : !!
2029 : !! FUNCTION
2030 : !! Generate a set of special k (or q) points which samples in a homogeneous way
2031 : !! the entire Brillouin zone of a simple lattice, face-centered cubic,
2032 : !! body-centered lattice and hexagonal lattice.
2033 : !! If kptrlatt is diagonal, the algorithm used here reduces to the usual
2034 : !! Monkhorst-Pack set of k points.
2035 : !!
2036 : !! INPUTS
2037 : !! brav = 1 or -1 -> simple lattice; 2 -> face-centered cubic;
2038 : !! 3 -> body-centered lattice; 4 -> hexagonal lattice (D6h)
2039 : !! downsampling(3) [optional, for brav=1 only]
2040 : !! Three integer numbers, describing the downsampling of the k grid
2041 : !! If present, in any case, only the first shiftk is taken into account
2042 : !! The absolute value of one number gives, for the corresponding k-coordinate, the factor of decrease of the sampling
2043 : !! If zero, only one point is used to sample along this direction
2044 : !! The sign has also a meaning :
2045 : !! - if three numbers are negative, perform a face-centered sampling
2046 : !! - if two numbers are negative, perform a body-centered sampling
2047 : !! - if one number is negative, perform a face-centered sampling for the two-dimensional lattice of the other directions
2048 : !! - if one number is zero and at least one number is negative, perform face-centered sampling for the non-zero directions.
2049 : !! iout = unit number for output
2050 : !! kptrlatt(3,3)=integer coordinates of the primitive vectors of the
2051 : !! lattice reciprocal to the k point lattice to be generated here
2052 : !! If diagonal, the three values are the Monkhorst-Pack usual values, in case of simple cubic.
2053 : !! mkpt = maximum number of k points
2054 : !! nshiftk= number of shift vectors in the repeated cell
2055 : !! option= Flag defining what will be printed of iout: 0 for k points, anything else for q points.
2056 : !! Also, for q points, if the Gamma point is present, place it first in the list.
2057 : !! shiftk(3,nshiftk) = vectors that will be used to determine the shifts from (0. 0. 0.).
2058 : !! prtout=write info to the output and log files.
2059 : !!
2060 : !! OUTPUT
2061 : !! nkpt = number of k points
2062 : !! spkpt(3,mkpt) = the nkpt first values contain the special k points
2063 : !! obtained by the Monkhorst & Pack method, in reduced coordinates.
2064 : !! These vectors have to be multiplied by the reciprocal basis vectors
2065 : !! gprimd(3,3) (in cartesian coordinates) to obtain the special k points
2066 : !! set in cartesian coordinates.
2067 : !!
2068 : !! NOTES
2069 : !! also allows for more than one vector in repeated cell.
2070 : !! this routine should be rewritten, to use the Wigner-Seitz cell,
2071 : !! and thus unify the different treatments.
2072 : !! References :
2073 : !! H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188 (1976) [[cite:Monkhorst1976]]
2074 : !! J.D. Pack and H.J. Monkhorst, Phys. Rev. B 16, 1748 (1977) [[cite:Pack1977]]
2075 : !! A.H. MacDonald, Phys. Rev. B 18, 5897 (1978) [[cite:MacDonald1978]]
2076 : !! R.A. Evarestov and V.P. Smirnov, Phys. Stat. Sol. (b) 119, 9 (1983) [[cite:Evarestov1983]]
2077 : !!
2078 : !! SOURCE
2079 :
2080 18071 : subroutine smpbz(brav,iout,kptrlatt,mkpt,nkpt,nshiftk,option,shiftk,spkpt,downsampling,prtout)
2081 :
2082 : !Arguments -------------------------------
2083 : !scalars
2084 : integer,intent(in) :: brav,iout,mkpt,nshiftk,option
2085 : integer,intent(out) :: nkpt
2086 : logical, optional, intent(in) :: prtout
2087 : !arrays
2088 : integer,intent(in) :: kptrlatt(3,3)
2089 : integer,optional,intent(in) :: downsampling(3)
2090 : real(dp),intent(in) :: shiftk(3,nshiftk)
2091 : real(dp),intent(out) :: spkpt(3,mkpt)
2092 :
2093 : !Local variables -------------------------
2094 : !scalars
2095 : integer,parameter :: prtvol=0
2096 : integer :: dividedown,ii,ikshft,jj,kk,nkpout,nkptlatt,nn,proddown
2097 : logical :: prtout_
2098 : real(dp) :: shift
2099 : character(len=500) :: msg
2100 : !arrays
2101 : integer :: ads(3),boundmax(3),boundmin(3),cds(3),coord(3),ngkpt(3)
2102 18071 : integer, allocatable :: found1(:,:),found2(:,:),found3(:,:)
2103 : real(dp) :: k1(3),k2(3),kcar(3),klatt(3,3),ktest(3),rlatt(3,3)
2104 : ! *********************************************************************
2105 :
2106 : !write(std_out,*)' smpbz : brav,iout,mkpt,nkpt,option=',brav,iout,mkpt,nkpt,option
2107 : !write(std_out,*)' smpbz : kptrlatt(:,:)=',kptrlatt(:,:)
2108 : !write(std_out,*)' smpbz : nshiftk=',nshiftk
2109 : !write(std_out,*)' smpbz : shiftk(:,:)=',shiftk(:,:)
2110 : !write(std_out,*)' smpbz : downsampling(:)=',downsampling(:)
2111 :
2112 18071 : prtout_ = .true.
2113 18071 : if (present(prtout)) then
2114 242 : prtout_ = prtout
2115 : end if
2116 :
2117 18071 : if (option/=0) then
2118 416 : if (prtout_) call wrtout(iout,' Homogeneous q point set in the B.Z. ')
2119 : end if
2120 :
2121 18071 : if(abs(brav)/=1)then
2122 : ! Only generate Monkhorst-Pack lattices
2123 : if(kptrlatt(1,2)/=0 .or. kptrlatt(2,1)/=0 .or. &
2124 : & kptrlatt(1,3)/=0 .or. kptrlatt(3,1)/=0 .or. &
2125 77 : & kptrlatt(2,3)/=0 .or. kptrlatt(3,2)/=0 ) then
2126 : write(msg, '(2a,a,3i0,a,a,3i4,a,a,3i4)' )&
2127 0 : & 'When abs(brav)/=1, kptrlatt must be diagonal, while it is',ch10,&
2128 0 : & 'kptrlatt(:,1)= ',kptrlatt(:,1),ch10,&
2129 0 : & 'kptrlatt(:,2)= ',kptrlatt(:,2),ch10,&
2130 0 : & 'kptrlatt(:,3)= ',kptrlatt(:,3)
2131 0 : ABI_BUG(msg)
2132 : end if
2133 :
2134 77 : ngkpt(1)=kptrlatt(1,1)
2135 77 : ngkpt(2)=kptrlatt(2,2)
2136 77 : ngkpt(3)=kptrlatt(3,3)
2137 : !
2138 77 : if( (ngkpt(1)<=0.or.ngkpt(2)<=0.or.ngkpt(3)<=0) .and. (ngkpt(1)/=0.or.ngkpt(2)/=0.or.ngkpt(3)/=0) ) then
2139 : write(msg, '(5a,i4,a,a,i0,a,a,i0,a,a)' )&
2140 0 : & 'All ngkpt (or ngqpt) must be strictly positive',ch10,&
2141 0 : & 'or all ngk(q)pt must be zero (for Gamma sampling), but :',ch10,&
2142 0 : & 'ngk(q)pt(1) = ',ngkpt(1),ch10,&
2143 0 : & 'ngk(q)pt(2) = ',ngkpt(2),ch10,&
2144 0 : & 'ngk(q)pt(3) = ',ngkpt(3),ch10,&
2145 0 : & 'Action: correct ngkpt or ngqpt in the input file.'
2146 0 : ABI_BUG(msg)
2147 : end if
2148 : end if
2149 :
2150 : !Just in case the user wants the grid downsampled to the Gamma point, checks that it is present, and possibly exits
2151 18071 : if(present(downsampling))then
2152 29064 : if(sum(abs(downsampling(:)))==0)then
2153 18 : do ikshft=1,nshiftk
2154 72 : if(sum(abs(shiftk(:,ikshft)))>tol12)cycle
2155 18 : nkpt=1
2156 72 : spkpt(:,1)=zero
2157 0 : return
2158 : end do
2159 : end if
2160 : end if
2161 :
2162 18053 : if(abs(brav)==1)then
2163 :
2164 : ! Compute the number of k points in the G-space unit cell
2165 : ! (will be multiplied by nshiftk later).
2166 : nkptlatt=kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3) &
2167 : & +kptrlatt(1,2)*kptrlatt(2,3)*kptrlatt(3,1) &
2168 : & +kptrlatt(1,3)*kptrlatt(2,1)*kptrlatt(3,2) &
2169 : & -kptrlatt(1,2)*kptrlatt(2,1)*kptrlatt(3,3) &
2170 : & -kptrlatt(1,3)*kptrlatt(2,2)*kptrlatt(3,1) &
2171 17976 : & -kptrlatt(1,1)*kptrlatt(2,3)*kptrlatt(3,2)
2172 :
2173 17976 : if(present(downsampling))then
2174 7248 : if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then
2175 150 : if(nshiftk>1)then
2176 : write(msg, '(a,3i4,2a,i4,4a)' )&
2177 0 : & 'Real downsampling is activated, with downsampling(1:3)=',downsampling(1:3),ch10,&
2178 0 : & 'However, nshiftk must be 1 in this case, while the input nshiftk=',nshiftk,ch10,&
2179 0 : & 'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,&
2180 0 : & 'or set nshiftk=1.'
2181 0 : ABI_ERROR(msg)
2182 : end if
2183 150 : proddown=downsampling(1)*downsampling(2)*downsampling(3)
2184 150 : if(proddown/=0)then
2185 150 : dividedown=abs(proddown)
2186 600 : if(minval(downsampling(:))<0)then ! If there is at least one negative number
2187 42 : dividedown=dividedown*2
2188 42 : if(proddown>0)dividedown=dividedown*2 ! If there are two negative numbers
2189 : end if
2190 : end if
2191 150 : if(mod(nkptlatt,dividedown)==0)then
2192 150 : nkptlatt=nkptlatt/dividedown
2193 : else
2194 : write(msg, '(a,3i4,2a,i4,4a)' )&
2195 0 : & 'The requested downsampling, with downsampling(1:3)=',downsampling(1:3),ch10,&
2196 0 : & 'is not compatible with kptrlatt=',ch10,&
2197 0 : & kptrlatt(:,:),ch10,&
2198 0 : & 'that gives nkptlatt=',nkptlatt,ch10,&
2199 0 : & 'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,&
2200 0 : & 'or modify your k-point grid and/or your downsampling in order for them to be compatible.'
2201 0 : ABI_ERROR(msg)
2202 : end if
2203 : end if
2204 : end if
2205 :
2206 : ! Simple Lattice
2207 : if (prtvol > 0) call wrtout(std_out,' Simple Lattice Grid ')
2208 17976 : if (mkpt<nkptlatt*nshiftk) then
2209 : write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
2210 0 : & 'The value of mkpt is not large enough. It should be',ch10,&
2211 0 : & 'at least',nkptlatt*nshiftk,',',ch10,&
2212 0 : & 'Action: set mkpt to that value in the main routine,',ch10,&
2213 0 : & 'and recompile the code.'
2214 0 : ABI_BUG(msg)
2215 : end if
2216 :
2217 : ! Build primitive vectors of the k lattice
2218 233688 : rlatt(:,:)=kptrlatt(:,:)
2219 17976 : call matr3inv(rlatt,klatt)
2220 :
2221 : ! write(std_out,*)' First primitive vector of the k lattice :',klatt(:,1)
2222 : ! write(std_out,*)' Second primitive vector of the k lattice :',klatt(:,2)
2223 : ! write(std_out,*)' Third primitive vector of the k lattice :',klatt(:,3)
2224 :
2225 : ! Now, klatt contains the three primitive vectors of the k lattice,
2226 : ! in reduced coordinates. One builds all k vectors that
2227 : ! are contained in the first Brillouin zone, with coordinates
2228 : ! in the interval [0,1[ . First generate boundaries of a big box.
2229 :
2230 71904 : do jj=1,3
2231 :
2232 : ! Mathematically, one has to find the coordinates of the corners of a
2233 : ! rectangular paralleliped with integer coordinates, that multiplies the klatt primitive cell and allows
2234 : ! it to incorporate completely the [0,1]^3 box. Then take the minimum and maximum
2235 : ! of these coordinates, and round them negatively and positively to the next integer.
2236 : ! This can be done easily using kptrlatt, considering each coordinate in turn
2237 : ! and boils down to enlarging the boundaries for jj by the value of kptrlatt(:,jj),
2238 : ! acting on boundmin or boundmax depending on the sign of kptrlatt(:,jj).
2239 : ! XG171020 The coding before 171020 was correct, despite being very simple.
2240 53928 : boundmin(jj)=0 ; boundmax(jj)=0
2241 215712 : do ii=1,3
2242 161784 : if(kptrlatt(ii,jj)<0)boundmin(jj)=boundmin(jj)+kptrlatt(ii,jj)
2243 215712 : if(kptrlatt(ii,jj)>0)boundmax(jj)=boundmax(jj)+kptrlatt(ii,jj)
2244 : end do
2245 :
2246 : ! To accommodate the shifts, boundmin and boundmax don't start from 0, but are enlarged by one
2247 : ! positively and/or negatively.
2248 : ! XG171020 Coding in v8.6.0 and before was not correct. This one is even simpler actually.
2249 161805 : boundmin(jj)=boundmin(jj)-ceiling(maxval(shiftk(jj,:))+tol14)
2250 179781 : boundmax(jj)=boundmax(jj)-floor(minval(shiftk(jj,:))-tol14)
2251 :
2252 : end do
2253 :
2254 17976 : if(present(downsampling))then
2255 28992 : ABI_MALLOC(found1,(boundmin(2):boundmax(2),boundmin(3):boundmax(3)))
2256 28992 : ABI_MALLOC(found2,(boundmin(1):boundmax(1),boundmin(3):boundmax(3)))
2257 28992 : ABI_MALLOC(found3,(boundmin(1):boundmax(1),boundmin(2):boundmax(2)))
2258 1983234 : found1=0 ; found2=0 ; found3=0
2259 : end if
2260 :
2261 17976 : nn=1
2262 164086 : do kk=boundmin(3),boundmax(3)
2263 146110 : coord(3)=kk
2264 1845684 : do jj=boundmin(2),boundmax(2)
2265 1681598 : coord(2)=jj
2266 29674069 : do ii=boundmin(1),boundmax(1)
2267 27846361 : coord(1)=ii
2268 :
2269 : ! Here, apply the downsampling : skip some of the trials
2270 27846361 : if(present(downsampling))then
2271 :
2272 9279498 : if(downsampling(1)==0 .and. found1(coord(2),coord(3))==1)cycle
2273 9279498 : if(downsampling(2)==0 .and. found2(coord(1),coord(3))==1)cycle
2274 9279498 : if(downsampling(3)==0 .and. found3(coord(1),coord(2))==1)cycle
2275 :
2276 37117992 : ads(:)=abs(downsampling(:))
2277 9279498 : if(ads(1)>0 .and. mod(coord(1),ads(1))/=0)cycle
2278 9230166 : if(ads(2)>0 .and. mod(coord(2),ads(2))/=0)cycle
2279 9208398 : if(ads(3)>0 .and. mod(coord(3),ads(2))/=0)cycle
2280 36795160 : cds(:)=coord(:)/ads(:)
2281 36795160 : if(minval(downsampling(:))<0)then ! If there is at least one negative number
2282 :
2283 1134 : if(downsampling(1)*downsampling(2)*downsampling(3)/=0)then ! If there is no zero number
2284 : ! Face-centered case
2285 1134 : if(downsampling(1)<0 .and. downsampling(2)<0 .and. downsampling(3)<0)then ! All three are negative
2286 1512 : if(mod(sum(cds(:)),2)/=0)cycle
2287 : ! One-face-centered case
2288 756 : else if(downsampling(1)*downsampling(2)*downsampling(3)<0)then ! Only one is negative
2289 378 : if(downsampling(1)<0 .and. mod(cds(2)+cds(3),2)/=0)cycle
2290 210 : if(downsampling(2)<0 .and. mod(cds(1)+cds(3),2)/=0)cycle
2291 210 : if(downsampling(3)<0 .and. mod(cds(1)+cds(2),2)/=0)cycle
2292 : ! Body-centered case ! What is left : two are negative
2293 : else
2294 : ! Either all are zero, or all are one, so skip when sum is 1 or 2.
2295 2646 : if(sum(mod(cds(:),2))==1 .or. sum(mod(cds(:),2))==2)cycle
2296 : end if
2297 : else
2298 0 : if(downsampling(1)==0 .and. mod(cds(2)+cds(3),2)/=0)cycle
2299 0 : if(downsampling(2)==0 .and. mod(cds(1)+cds(3),2)/=0)cycle
2300 0 : if(downsampling(3)==0 .and. mod(cds(1)+cds(2),2)/=0)cycle
2301 : end if
2302 : end if
2303 : end if
2304 :
2305 57212392 : do ikshft=1,nshiftk
2306 :
2307 : ! Only the first shiftk is taken into account if downsampling
2308 : ! if(.false.)then
2309 27765743 : if(present(downsampling))then
2310 9198188 : if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then
2311 6996 : if(ikshft>1)cycle
2312 : end if
2313 : end if
2314 :
2315 : ! Coordinates of the trial k point with respect to the k primitive lattice
2316 27765743 : k1(1)=ii+shiftk(1,ikshft)
2317 27765743 : k1(2)=jj+shiftk(2,ikshft)
2318 27765743 : k1(3)=kk+shiftk(3,ikshft)
2319 : ! Reduced coordinates of the trial k point
2320 111062972 : k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3)
2321 : ! Eliminate the point if outside [0,1[
2322 27765743 : if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle
2323 16527426 : if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle
2324 10898485 : if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle
2325 : ! Wrap the trial values in the interval ]-1/2,1/2] .
2326 7075194 : call wrap2_pmhalf(k2(1),k1(1),shift)
2327 7075194 : call wrap2_pmhalf(k2(2),k1(2),shift)
2328 7075194 : call wrap2_pmhalf(k2(3),k1(3),shift)
2329 28300776 : spkpt(:,nn)=k1(:)
2330 7075194 : nn=nn+1
2331 :
2332 34921555 : if(present(downsampling))then
2333 1482698 : found1(coord(2),coord(3))=1
2334 1482698 : found2(coord(1),coord(3))=1
2335 1482698 : found3(coord(1),coord(2))=1
2336 : end if
2337 :
2338 : end do
2339 : end do
2340 : end do
2341 : end do
2342 17976 : nkpt=nn-1
2343 :
2344 17976 : if(present(downsampling))then
2345 7248 : ABI_FREE(found1)
2346 7248 : ABI_FREE(found2)
2347 7248 : ABI_FREE(found3)
2348 : end if
2349 :
2350 17976 : if(nkpt/=nkptlatt*nshiftk)then
2351 : write(msg, '(a,i0,3a,i0,a)' )&
2352 0 : 'The number of k points ',nkpt,' is not equal to',ch10,&
2353 0 : 'nkptlatt*nshiftk which is ',nkptlatt*nshiftk,'.'
2354 0 : ABI_BUG(msg)
2355 : end if
2356 :
2357 77 : else if(brav==2)then
2358 :
2359 : ! Face-Centered Lattice
2360 : if (prtvol > 0) call wrtout(std_out,' Face-Centered Lattice Grid ')
2361 49 : if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2) then
2362 : write(msg, '(a,a,a,i0,a,a,a,a,a)' )&
2363 0 : & 'The value of mkpt is not large enough. It should be',ch10,&
2364 0 : & 'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,',',ch10,&
2365 0 : & 'Action: set mkpt to that value in the main routine,',ch10,&
2366 0 : & 'and recompile the code.'
2367 0 : ABI_BUG(msg)
2368 : end if
2369 49 : nn=1
2370 49 : if (ngkpt(1)/=ngkpt(2).or.ngkpt(1)/=ngkpt(3)) then
2371 : write(msg, '(4a,3(a,i0,a),a)' )&
2372 0 : & 'For face-centered lattices, the numbers ngqpt(1:3)',ch10,&
2373 0 : & 'must be equal, while they are :',ch10,&
2374 0 : & 'ngqpt(1) = ',ngkpt(1),ch10,&
2375 0 : & 'ngqpt(2) = ',ngkpt(2),ch10,&
2376 0 : & 'ngqpt(3) = ',ngkpt(3),ch10,&
2377 0 : & 'Action: modify ngqpt(1:3) in the input file.'
2378 0 : ABI_BUG(msg)
2379 : end if
2380 49 : if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2)) then
2381 : write(msg, '(4a,3(a,i0,a),a)' )&
2382 0 : & 'For face-centered lattices, the numbers ngqpt(1:3)*nshiftk',ch10,&
2383 0 : & 'must be even, while they are :',ch10,&
2384 0 : & 'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,&
2385 0 : & 'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,&
2386 0 : & 'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,&
2387 0 : & 'Action: modify ngqpt(1:3)*nshiftk in the input file.'
2388 0 : ABI_ERROR(msg)
2389 : end if
2390 49 : if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then
2391 0 : spkpt(1,1)=0.0_dp
2392 0 : spkpt(2,1)=0.0_dp
2393 0 : spkpt(3,1)=0.0_dp
2394 0 : nkpt=1
2395 : else
2396 390 : do kk=1,ngkpt(3)
2397 3223 : do jj=1,ngkpt(2)
2398 31751 : do ii=1,ngkpt(1)
2399 60204 : do ikshft=1,nshiftk
2400 28794 : k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1)
2401 28794 : k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2)
2402 28794 : k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3)
2403 : ! Wrap the trial values in the interval ]-1/2,1/2] .
2404 28794 : call wrap2_pmhalf(k1(1),k2(1),shift)
2405 28794 : call wrap2_pmhalf(k1(2),k2(2),shift)
2406 28794 : call wrap2_pmhalf(k1(3),k2(3),shift)
2407 : ! Test whether it is inside the FCC BZ.
2408 28794 : ktest(1)=2*k2(1)-1.0d-10
2409 28794 : ktest(2)=2*k2(2)-2.0d-10
2410 28794 : ktest(3)=2*k2(3)-5.0d-10
2411 57371 : if (abs(ktest(1))+abs(ktest(2))+abs(ktest(3))<1.5_dp) then
2412 14397 : kcar(1)=ktest(1)+1.0d-10
2413 14397 : kcar(2)=ktest(2)+2.0d-10
2414 14397 : kcar(3)=ktest(3)+5.0d-10
2415 14397 : spkpt(1,nn)=0.5_dp*kcar(2)+0.5_dp*kcar(3)
2416 14397 : spkpt(2,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(3)
2417 14397 : spkpt(3,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(2)
2418 14397 : nn=nn+1
2419 : end if
2420 : end do
2421 : end do
2422 : end do
2423 : end do
2424 49 : nkpt=nn-1
2425 49 : if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2)then
2426 : write(msg, '(a,i8,a,a,a,i8,a)' )&
2427 0 : & 'The number of k points ',nkpt,' is not equal to',ch10,&
2428 0 : & '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2 which is',&
2429 0 : & (ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,'.'
2430 0 : ABI_BUG(msg)
2431 : end if
2432 : end if
2433 :
2434 28 : else if(brav==3)then
2435 :
2436 : ! Body-Centered Lattice (not mandatory cubic !)
2437 : if (prtvol > 0) call wrtout(std_out,' Body-Centered Lattice Grid ')
2438 13 : if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/4) then
2439 : write(msg, '(a,a,a,i8,a,a,a,a,a)' )&
2440 0 : & 'The value of mkpt is not large enough. It should be',ch10,&
2441 0 : & 'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,',',ch10,&
2442 0 : & 'Action: set mkpt to that value in the main routine,',ch10,&
2443 0 : & 'and recompile the code.'
2444 0 : ABI_BUG(msg)
2445 : end if
2446 13 : nn=1
2447 : if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2) .or.&
2448 13 : & (ngkpt(2)*nshiftk)/=(((ngkpt(2)*nshiftk)/2)*2) .or.&
2449 : & (ngkpt(3)*nshiftk)/=(((ngkpt(3)*nshiftk)/2)*2) ) then
2450 : write(msg, '(4a,3(a,i6,a),a)' )&
2451 0 : & 'For body-centered lattices, the numbers ngqpt(1:3)',ch10,&
2452 0 : & 'must be even, while they are :',ch10,&
2453 0 : & 'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,&
2454 0 : & 'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,&
2455 0 : & 'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,&
2456 0 : & 'Action: modify ngqpt(1:3) in the input file.'
2457 0 : ABI_ERROR(msg)
2458 : end if
2459 13 : if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then
2460 0 : spkpt(1,1)=0.0_dp
2461 0 : spkpt(2,1)=0.0_dp
2462 0 : spkpt(3,1)=0.0_dp
2463 0 : nkpt=1
2464 : else
2465 117 : do kk=1,ngkpt(3)
2466 949 : do jj=1,ngkpt(2)
2467 7592 : do ii=1,ngkpt(1)
2468 14144 : do ikshft=1,nshiftk
2469 6656 : k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1)
2470 6656 : k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2)
2471 6656 : k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3)
2472 : ! Wrap the trial values in the interval ]-1/2,1/2] .
2473 6656 : call wrap2_pmhalf(k1(1),k2(1),shift)
2474 6656 : call wrap2_pmhalf(k1(2),k2(2),shift)
2475 6656 : call wrap2_pmhalf(k1(3),k2(3),shift)
2476 : ! Test whether it is inside the BCC BZ.
2477 6656 : ktest(1)=2*k2(1)-1.0d-10
2478 6656 : ktest(2)=2*k2(2)-2.0d-10
2479 6656 : ktest(3)=2*k2(3)-5.0d-10
2480 13312 : if (abs(ktest(1))+abs(ktest(2))<1._dp) then
2481 3328 : if (abs(ktest(1))+abs(ktest(3))<1._dp) then
2482 2184 : if (abs(ktest(2))+abs(ktest(3))<1._dp) then
2483 1664 : kcar(1)=ktest(1)+1.0d-10
2484 1664 : kcar(2)=ktest(2)+2.0d-10
2485 1664 : kcar(3)=ktest(3)+5.0d-10
2486 1664 : spkpt(1,nn)=-0.5*kcar(1)+0.5*kcar(2)+0.5*kcar(3)
2487 1664 : spkpt(2,nn)=0.5*kcar(1)-0.5*kcar(2)+0.5*kcar(3)
2488 1664 : spkpt(3,nn)=0.5*kcar(1)+0.5*kcar(2)-0.5*kcar(3)
2489 1664 : nn=nn+1
2490 : end if
2491 : end if
2492 : end if
2493 : end do
2494 : end do
2495 : end do
2496 : end do
2497 13 : nkpt=nn-1
2498 13 : if(nkpt==0)then
2499 : write(msg, '(3a)' )&
2500 0 : & 'BCC lattice, input ngqpt=0, so no kpt is generated.',ch10,&
2501 0 : & 'Action: modify ngqpt(1:3) in the input file.'
2502 0 : ABI_ERROR(msg)
2503 : end if
2504 13 : if(nkpt/=(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4)then
2505 : write(msg, '(a,i0,3a,i0,a)' )&
2506 0 : & 'The number of k points ',nkpt,' is not equal to',ch10,&
2507 0 : & '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4 which is',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,'.'
2508 0 : ABI_BUG(msg)
2509 : end if
2510 : end if
2511 :
2512 15 : else if(brav==4)then
2513 :
2514 : ! Hexagonal Lattice (D6h)
2515 : if (prtvol > 0) call wrtout(std_out,' Hexagonal Lattice Grid ')
2516 15 : if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)) then
2517 : write(msg, '(a,a,a,i0,a,a,a,a,a)' )&
2518 0 : & 'The value of mkpt is not large enough. It should be',ch10,&
2519 0 : & 'at least',ngkpt(1)*ngkpt(2)*ngkpt(3),',',ch10,&
2520 0 : & 'Action: set mkpt to that value in the main routine,',ch10,&
2521 0 : & 'and recompile the code.'
2522 0 : ABI_BUG(msg)
2523 : end if
2524 15 : nn=1
2525 15 : if (ngkpt(1)/=ngkpt(2)) then
2526 : write(msg, '(4a,2(a,i0,a),a)' )&
2527 0 : & 'For hexagonal lattices, the numbers ngqpt(1:2)',ch10,&
2528 0 : & 'must be equal, while they are:',ch10,&
2529 0 : & 'ngqpt(1) = ',ngkpt(1),ch10,&
2530 0 : & 'ngqpt(2) = ',ngkpt(2),ch10,&
2531 0 : & 'Action: modify ngqpt(1:3) in the input file.'
2532 0 : ABI_ERROR(msg)
2533 : end if
2534 15 : if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then
2535 : write(msg, '(3a)' )&
2536 0 : & 'For hexagonal lattices, ngqpt(1:3)=0 is not permitted',ch10,&
2537 0 : & 'Action: modify ngqpt(1:3) in the input file.'
2538 0 : ABI_ERROR(msg)
2539 : else
2540 120 : do kk=1,ngkpt(3)
2541 1161 : do jj=1,ngkpt(2)
2542 14459 : do ii=1,ngkpt(1)
2543 27667 : do ikshft=1,nshiftk
2544 13313 : k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1)
2545 13313 : k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2)
2546 13313 : k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3)
2547 : ! Wrap the trial values in the interval ]-1/2,1/2] .
2548 13313 : call wrap2_pmhalf(k1(1),k2(1),shift)
2549 13313 : call wrap2_pmhalf(k1(2),k2(2),shift)
2550 13313 : call wrap2_pmhalf(k1(3),k2(3),shift)
2551 53252 : spkpt(:,nn)=k2(:)
2552 26626 : nn=nn+1
2553 : end do
2554 : end do
2555 : end do
2556 : end do
2557 15 : nkpt=nn-1
2558 15 : if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)then
2559 : write(msg, '(a,i0,3a,i0,a)' )&
2560 0 : & 'The number of k points ',nkpt,' is not equal to',ch10,&
2561 0 : & 'ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk which is',ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk,'.'
2562 0 : ABI_BUG(msg)
2563 : end if
2564 : end if
2565 :
2566 : else
2567 :
2568 : write(msg, '(a,i0,a,a,a)' )&
2569 0 : & 'The calling routine asks brav= ',brav,'.',ch10,&
2570 0 : & 'but only brav=1 or -1,2,3 or 4 are allowed.'
2571 0 : ABI_BUG(msg)
2572 : end if
2573 :
2574 18053 : if (option/=0) then
2575 : ! Put the Gamma point first
2576 416 : if(nkpt>1)then
2577 45486 : do ii=1,nkpt
2578 181662 : if(sum(abs(spkpt(:,ii)))<tol8)then
2579 1192 : spkpt(:,ii)=spkpt(:,1)
2580 1192 : spkpt(:,1)=zero
2581 : exit
2582 : end if
2583 : end do
2584 : end if
2585 :
2586 416 : if (prtout_) then
2587 370 : write(msg,'(a,i8)')' Grid q points : ',nkpt
2588 370 : call wrtout(iout,msg)
2589 370 : nkpout=nkpt
2590 370 : if(nkpt>80)then
2591 77 : call wrtout(iout,' greater than 80, so only write 20 of them ')
2592 77 : nkpout=20
2593 : end if
2594 8991 : do ii=1,nkpout
2595 8621 : write(msg, '(1x,i2,a2,3es16.8)' )ii,') ',spkpt(1,ii),spkpt(2,ii),spkpt(3,ii)
2596 8991 : call wrtout(iout,msg)
2597 : end do
2598 : end if
2599 : end if
2600 :
2601 18053 : end subroutine smpbz
2602 : !!***
2603 :
2604 : !!****f* m_kpts/testkgrid
2605 : !! NAME
2606 : !! testkgrid
2607 : !!
2608 : !! FUNCTION
2609 : !! Test different grids of k points. The algorithm used is based on the idea of testing different
2610 : !! one-dimensional sets of possible k point grids. It is not exhaustive (other families could be included),
2611 : !! but should do a respectable job in all cases. The Monkhorst-Pack set of grids (defined with respect to
2612 : !! symmetry axes, and not primitive axes) is always tested.
2613 : !!
2614 : !! INPUTS
2615 : !! bravais(11): bravais(1)=iholohedry
2616 : !! bravais(2)=center
2617 : !! bravais(3:11)=coordinates of rprim in the axes of the conventional bravais lattice (*2 if center/=0)
2618 : !! iout=unit number for echoed output
2619 : !! msym=default maximal number of symmetries
2620 : !! nsym=number of symmetries
2621 : !! prtkpt=if non-zero, will write the characteristics of k grids, then stop
2622 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
2623 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
2624 : !! symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
2625 : !! vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
2626 : !!
2627 : !! OUTPUT
2628 : !! kptrlatt(3,3)=k-point lattice specification
2629 : !! nshiftk=number of k-point shifts in shiftk (always 1 from this routine)
2630 : !! shiftk(3,MAX_NSHIFTK)=shift vectors for k point generation
2631 : !!
2632 : !! SIDE EFFECTS
2633 : !! kptrlen=length of the smallest real space supercell vector associated with the lattice of k points.
2634 : !!
2635 : !! NOTES
2636 : !! Note that nkpt can be computed by calling this routine with input value nkpt=0
2637 : !! Note that kptopt is always =1 in this routine.
2638 : !!
2639 : !! SOURCE
2640 :
2641 197 : subroutine testkgrid(bravais,iout,kptrlatt,kptrlen,msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum)
2642 :
2643 : !Arguments ------------------------------------
2644 : !scalars
2645 : integer,intent(in) :: iout,msym,nsym,prtkpt
2646 : integer,intent(out) :: nshiftk
2647 : real(dp),intent(inout) :: kptrlen
2648 : !arrays
2649 : integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3)
2650 : integer,intent(out) :: kptrlatt(3,3)
2651 : real(dp),intent(in) :: rprimd(3,3)
2652 : real(dp),intent(inout) :: shiftk(3,MAX_NSHIFTK) !vz_i
2653 :
2654 : !Local variables-------------------------------
2655 : !scalars
2656 : integer,parameter :: kptopt=1,mkpt_list=100000
2657 : integer :: ang90,center,dirvacuum,equal,igrid,igrid_current,iholohedry,ii,init_mult,iscale,iscf
2658 : integer :: iset,mult1,mult2,mult3,ndims,nkpt,nkpt_current,nkpt_trial,nset
2659 : real(dp) :: buffer_scale,determinant,fact,factor,kptrlen_current,kptrlen_max,kptrlen_target
2660 : real(dp) :: kptrlen_trial,length1,length2,length3,length_axis1,length_axis2
2661 : real(dp) :: length_axis3,merit_factor,mult1h,mult2h,mult3h,reduceda,reducedb
2662 : real(dp) :: sca,scb,scc,surface,ucvol
2663 : character(len=500) :: msg
2664 : !arrays
2665 : integer :: kptrlatt_current(3,3),kptrlatt_trial(3,3)
2666 197 : integer,allocatable :: grid_list(:)
2667 : real(dp) :: axes(3,3),gmet(3,3),gprimd(3,3),matrix1(3,3),matrix2(3,3)
2668 : real(dp) :: metmin(3,3),minim(3,3),r2d(3,3),rmet(3,3),rsuper(3,3)
2669 : real(dp) :: shiftk_current(3,MAX_NSHIFTK),shiftk_trial(3,MAX_NSHIFTK)
2670 197 : real(dp),allocatable :: kpt(:,:),kptrlen_list(:),wtk(:)
2671 : ! *************************************************************************
2672 :
2673 197 : kptrlen_target=kptrlen
2674 :
2675 : !The vacuum array must be made of 0 or 1
2676 788 : do ii=1,3
2677 788 : if(vacuum(ii)/=0 .and. vacuum(ii)/=1)then
2678 : write(msg,'(a,a,a,i1,a,i3,a,a)')&
2679 0 : & 'The values of vacuum must be 0 or 1.',ch10,&
2680 0 : & 'However, the input vacuum(',ii,') is',vacuum(ii),ch10,&
2681 0 : & 'Action: correct vacuum in your input file.'
2682 0 : ABI_ERROR(msg)
2683 : end if
2684 : end do
2685 :
2686 : !Specific preparation for 2-dimensional system
2687 788 : if(sum(vacuum(:))==1)then
2688 :
2689 : ! Make the non-active vector orthogonal to the active vectors,
2690 : ! and take it along the z direction
2691 8 : if(vacuum(1)==1)then
2692 2 : r2d(1,3)=rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3)
2693 2 : r2d(2,3)=rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3)
2694 2 : r2d(3,3)=rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3)
2695 8 : r2d(:,1)=rprimd(:,2)
2696 8 : r2d(:,2)=rprimd(:,3)
2697 : dirvacuum=1
2698 6 : else if(vacuum(2)==1)then
2699 0 : r2d(1,3)=rprimd(2,3)*rprimd(3,1)-rprimd(3,3)*rprimd(2,1)
2700 0 : r2d(2,3)=rprimd(3,3)*rprimd(1,1)-rprimd(1,3)*rprimd(3,1)
2701 0 : r2d(3,3)=rprimd(1,3)*rprimd(2,1)-rprimd(2,3)*rprimd(1,1)
2702 0 : r2d(:,1)=rprimd(:,3)
2703 0 : r2d(:,2)=rprimd(:,1)
2704 : dirvacuum=2
2705 6 : else if(vacuum(3)==1)then
2706 6 : r2d(1,3)=rprimd(2,1)*rprimd(3,2)-rprimd(3,1)*rprimd(2,2)
2707 6 : r2d(2,3)=rprimd(3,1)*rprimd(1,2)-rprimd(1,1)*rprimd(3,2)
2708 6 : r2d(3,3)=rprimd(1,1)*rprimd(2,2)-rprimd(2,1)*rprimd(1,2)
2709 24 : r2d(:,1)=rprimd(:,1)
2710 24 : r2d(:,2)=rprimd(:,2)
2711 : dirvacuum=3
2712 : end if
2713 32 : surface=sqrt(sum(r2d(:,3)**2))
2714 : ! Identify the 2-D Bravais lattice
2715 : ! DEBUG
2716 : ! write(std_out,*)' r2d=',r2d(:,:)
2717 : ! ENDDEBUG
2718 8 : call metric(gmet,gprimd,-1,rmet,r2d,ucvol)
2719 8 : call smallprim(metmin,minim,r2d)
2720 : ! DEBUG
2721 : ! write(std_out,*)' minim=',minim(:,:)
2722 : ! ENDDEBUG
2723 8 : ang90=0 ; equal=0 ; center=0
2724 8 : axes(:,:)=minim(:,:)
2725 8 : if(abs(metmin(1,2))<tol8)ang90=1
2726 8 : if(abs(metmin(1,1)-metmin(2,2))<tol8)equal=1
2727 8 : if(ang90==1)then
2728 : if(equal==1)iholohedry=4
2729 5 : if(equal==0)iholohedry=2
2730 3 : else if(equal==1)then
2731 2 : reduceda=metmin(1,2)/metmin(1,1)
2732 2 : if(abs(reduceda+0.5_dp)<tol8)then
2733 : iholohedry=3
2734 1 : else if(abs(reduceda-0.5_dp)<tol8)then
2735 0 : iholohedry=3
2736 : ! Use conventional axes
2737 0 : axes(:,2)=minim(:,2)-minim(:,1)
2738 : else
2739 4 : iholohedry=2 ; center=1
2740 4 : axes(:,1)=minim(:,1)+minim(:,2)
2741 4 : axes(:,2)=minim(:,2)-minim(:,1)
2742 : end if
2743 : else
2744 1 : reduceda=metmin(1,2)/metmin(1,1)
2745 1 : reducedb=metmin(1,2)/metmin(2,2)
2746 1 : if(abs(reduceda+0.5_dp)<tol8)then
2747 0 : iholohedry=2 ; center=1
2748 0 : axes(:,2)=2.0_dp*minim(:,2)+minim(:,1)
2749 1 : else if(abs(reduceda-0.5_dp)<tol8)then
2750 0 : iholohedry=2 ; center=1
2751 0 : axes(:,2)=2.0_dp*minim(:,2)-minim(:,1)
2752 1 : else if(abs(reducedb+0.5_dp)<tol8)then
2753 0 : iholohedry=2 ; center=1
2754 0 : axes(:,1)=2.0_dp*minim(:,1)+minim(:,2)
2755 1 : else if(abs(reducedb-0.5_dp)<tol8)then
2756 0 : iholohedry=2 ; center=1
2757 0 : axes(:,1)=2.0_dp*minim(:,1)-minim(:,2)
2758 : else
2759 : iholohedry=1
2760 : end if
2761 : end if
2762 : ! Make sure that axes form a right-handed coordinate system
2763 : determinant=axes(1,1)*axes(2,2)*axes(3,3) &
2764 : & +axes(1,2)*axes(2,3)*axes(3,1) &
2765 : & +axes(1,3)*axes(3,2)*axes(2,1) &
2766 : & -axes(1,1)*axes(3,2)*axes(2,3) &
2767 : & -axes(1,3)*axes(2,2)*axes(3,1) &
2768 8 : & -axes(1,2)*axes(2,1)*axes(3,3)
2769 8 : if(determinant<zero)then
2770 0 : axes(:,1)=-axes(:,1)
2771 : end if
2772 : ! Prefer symmetry axes on the same side as the primitive axes
2773 32 : sca=DOT_PRODUCT(axes(:,1), r2d(:,1))
2774 32 : scb=DOT_PRODUCT(axes(:,2), r2d(:,2))
2775 : scc=DOT_PRODUCT(axes(:,3), rprimd(:,dirvacuum))
2776 8 : if(sca<-tol8 .and. scb<-tol8)then
2777 0 : axes(:,1)=-axes(:,1) ; sca=-sca
2778 8 : axes(:,2)=-axes(:,2) ; scb=-scb
2779 : end if
2780 : ! Doing this might change the angle between vectors, so that
2781 : ! the cell is not conventional anymore
2782 : ! if(sca<-tol8 .and. scc<-tol8)then
2783 : ! axes(:,1)=-axes(:,1) ; sca=-sca
2784 : ! axes(:,3)=-axes(:,3) ; scc=-scc
2785 : ! end if
2786 : ! if(scb<-tol8 .and. scc<-tol8)then
2787 : ! axes(:,2)=-axes(:,2) ; scb=-scb
2788 : ! axes(:,3)=-axes(:,3) ; scc=-scc
2789 : ! end if
2790 32 : length_axis1=NORM2(axes(:,1))
2791 32 : length_axis2=NORM2(axes(:,2))
2792 :
2793 : ! DEBUG
2794 : ! write(std_out,*)' testkgrid: iholohedry, center =',iholohedry,center
2795 : ! write(std_out,*)' testkgrid: axis 1=',axes(:,1)
2796 : ! write(std_out,*)' testkgrid: axis 2=',axes(:,2)
2797 : ! write(std_out,*)' testkgrid: axis 3=',axes(:,3)
2798 : ! write(std_out,*)' testkgrid: length_axis=',length_axis1,length_axis2
2799 : ! ENDDEBUG
2800 :
2801 : ! End special treatment of 2-D case
2802 : end if
2803 :
2804 : !3-dimensional system
2805 788 : if(sum(vacuum(:))==0)then
2806 83 : iholohedry=bravais(1)
2807 83 : center=bravais(2)
2808 83 : fact=1.0_dp
2809 83 : if(center/=0)fact=0.5_dp
2810 332 : matrix1(:,1)=bravais(3:5)*fact
2811 332 : matrix1(:,2)=bravais(6:8)*fact
2812 332 : matrix1(:,3)=bravais(9:11)*fact
2813 83 : call matr3inv(matrix1,matrix2)
2814 332 : do ii=1,3
2815 1079 : axes(:,ii)=rprimd(:,1)*matrix2(ii,1)+rprimd(:,2)*matrix2(ii,2)+rprimd(:,3)*matrix2(ii,3)
2816 : end do
2817 332 : length_axis1 = NORM2(axes(:,1))
2818 332 : length_axis2 = NORM2(axes(:,2))
2819 332 : length_axis3 = NORM2(axes(:,3))
2820 : ! DEBUG
2821 : ! write(std_out,*)' testkgrid: axes=',axes(:,:)
2822 : ! write(std_out,*)' length_axis=',length_axis1,length_axis2,length_axis3
2823 : ! ENDDEBUG
2824 : end if
2825 :
2826 : !This routine examine only primitive k lattices.
2827 197 : nshiftk=1
2828 :
2829 : !If prtkpt/=0, will examine more grids than strictly needed
2830 197 : buffer_scale=one
2831 197 : if(prtkpt/=0)buffer_scale=two
2832 :
2833 197 : if(prtkpt/=0)then
2834 13 : write(msg,'(a,a,a,a,a,a,a,a)' )ch10,&
2835 13 : ' testkgrid : will perform the analysis of a series of k-grids.',ch10,&
2836 13 : ' Note that kptopt=1 in this analysis, irrespective of its input value.',ch10,ch10,&
2837 26 : ' Grid# kptrlatt shiftk kptrlen nkpt iset',ch10
2838 13 : call wrtout(std_out,msg)
2839 13 : call wrtout(iout,msg)
2840 13 : ABI_MALLOC(grid_list,(mkpt_list))
2841 13 : ABI_MALLOC(kptrlen_list,(mkpt_list))
2842 1300013 : grid_list(:)=0
2843 1300013 : kptrlen_list(:)=0.0_dp
2844 : end if
2845 :
2846 788 : if(sum(vacuum(:))==3)then
2847 :
2848 102 : kptrlatt(:,:)=0
2849 102 : kptrlatt(1,1)=1
2850 102 : kptrlatt(2,2)=1
2851 102 : kptrlatt(3,3)=1
2852 408 : shiftk(:,1)=0.0_dp
2853 102 : kptrlen=1000.0_dp
2854 102 : nkpt_current=1
2855 102 : igrid_current=1
2856 :
2857 102 : if(prtkpt/=0)then
2858 : write(msg,&
2859 : & '(a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )&
2860 0 : & ' 1 ',kptrlatt(:,1),' ',shiftk(1,1),' ',kptrlen,1,1,ch10,&
2861 0 : & ' ',kptrlatt(:,2),' ',shiftk(2,1),ch10,&
2862 0 : & ' ',kptrlatt(:,3),' ',shiftk(3,1),ch10
2863 0 : call wrtout(std_out,msg)
2864 0 : call wrtout(iout,msg)
2865 : ! The unit cell volume is fake
2866 0 : ucvol=kptrlen**3
2867 : end if
2868 :
2869 : else
2870 :
2871 95 : nkpt=0 ; nkpt_current=0 ; iscf=1 ; iset=1
2872 95 : kptrlen_current=0.0_dp
2873 95 : mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1
2874 95 : ABI_MALLOC(kpt,(3,nkpt))
2875 95 : ABI_MALLOC(wtk,(nkpt))
2876 95 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2877 :
2878 : ! Loop on different grids, the upper limit is only to avoid an infinite loop
2879 2054 : do igrid=1,1000
2880 :
2881 2054 : kptrlatt_trial(:,:)=0
2882 2054 : kptrlatt_trial(1,1)=1
2883 2054 : kptrlatt_trial(2,2)=1
2884 2054 : kptrlatt_trial(3,3)=1
2885 8216 : shiftk_trial(:,1)=0.0_dp
2886 :
2887 : ! 1-dimensional system
2888 8216 : if(sum(vacuum(:))==2)then
2889 6 : if(vacuum(1)==0)then
2890 6 : kptrlatt_trial(1,1)=2*igrid ; shiftk_trial(1,1)=0.5_dp
2891 0 : else if(vacuum(2)==0)then
2892 0 : kptrlatt_trial(2,2)=2*igrid ; shiftk_trial(2,1)=0.5_dp
2893 0 : else if(vacuum(3)==0)then
2894 0 : kptrlatt_trial(3,3)=2*igrid ; shiftk_trial(3,1)=0.5_dp
2895 : end if
2896 : end if
2897 :
2898 : ! 2-dimensional system
2899 8216 : if(sum(vacuum(:))==1)then
2900 :
2901 : ! Treat hexagonal holohedries separately
2902 342 : if(iholohedry==3)then
2903 :
2904 : ! write(std_out,*)' testkgrid: 2D, hexagonal'
2905 :
2906 44 : mult1=mult1+1
2907 44 : nset=4
2908 44 : if(iset==1)then
2909 56 : rsuper(:,1)=axes(:,1)*mult1
2910 56 : rsuper(:,2)=axes(:,2)*mult1
2911 56 : shiftk_trial(:,1)=0.0_dp
2912 30 : else if(iset==2)then
2913 32 : rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1
2914 32 : rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
2915 8 : shiftk_trial(1,1)=1.0_dp/3.0_dp
2916 8 : shiftk_trial(2,1)=1.0_dp/3.0_dp
2917 22 : else if(iset==3)then
2918 32 : rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1
2919 32 : rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
2920 32 : shiftk_trial(:,1)=0.0_dp
2921 14 : else if(iset==4)then
2922 56 : rsuper(:,1)=axes(:,1)*mult1
2923 56 : rsuper(:,2)=axes(:,2)*mult1
2924 14 : shiftk_trial(1,1)=0.5_dp
2925 14 : shiftk_trial(2,1)=0.5_dp
2926 : end if
2927 :
2928 : else
2929 : ! Now treat all other holohedries
2930 298 : length1=length_axis1*mult1
2931 298 : length2=length_axis2*mult2
2932 : ! write(std_out,*)' testkgrid: (2d) length=',length1,length2
2933 298 : if(abs(length1-length2)<tol8)then
2934 103 : mult1=mult1+1
2935 103 : mult2=mult2+1
2936 195 : else if(length1>length2)then
2937 85 : mult2=mult2+1
2938 110 : else if(length2>length1)then
2939 110 : mult1=mult1+1
2940 : end if
2941 298 : nset=4
2942 : ! iset==5 and 6 are allowed only for centered lattice
2943 298 : if(center==1)nset=6
2944 298 : if(iset==1 .or. iset==2)then
2945 640 : rsuper(:,1)=axes(:,1)*mult1
2946 640 : rsuper(:,2)=axes(:,2)*mult2
2947 138 : else if(iset==3 .or. iset==4)then
2948 460 : rsuper(:,1)=axes(:,1)*mult1-axes(:,2)*mult2
2949 460 : rsuper(:,2)=axes(:,1)*mult1+axes(:,2)*mult2
2950 23 : else if(iset==5 .or. iset==6)then
2951 92 : rsuper(:,1)=axes(:,1)*(mult1-0.5_dp)-axes(:,2)*(mult2-0.5_dp)
2952 92 : rsuper(:,2)=axes(:,1)*(mult1-0.5_dp)+axes(:,2)*(mult2-0.5_dp)
2953 : end if
2954 : ! This was the easiest way to code all even mult1 and mult2 pairs:
2955 : ! make separate series for this possibility.
2956 298 : if(iset==2 .or. iset==4 .or. iset==6)then
2957 404 : rsuper(:,1)=2.0_dp*rsuper(:,1)
2958 404 : rsuper(:,2)=2.0_dp*rsuper(:,2)
2959 : end if
2960 298 : shiftk_trial(1,1)=0.5_dp
2961 298 : shiftk_trial(2,1)=0.5_dp
2962 :
2963 : end if
2964 :
2965 : ! Put back the inactive direction
2966 342 : if(dirvacuum==1)then
2967 56 : rsuper(:,3)=rsuper(:,1)
2968 14 : shiftk_trial(3,1)=shiftk_trial(1,1)
2969 56 : rsuper(:,1)=rprimd(:,1)
2970 14 : shiftk_trial(1,1)=0.0_dp
2971 328 : else if(dirvacuum==2)then
2972 0 : rsuper(:,3)=rsuper(:,1)
2973 0 : shiftk_trial(3,1)=shiftk_trial(1,1)
2974 0 : rsuper(:,1)=rsuper(:,2)
2975 0 : shiftk_trial(1,1)=shiftk_trial(2,1)
2976 0 : rsuper(:,2)=rprimd(:,2)
2977 0 : shiftk_trial(2,1)=0.0_dp
2978 : else if(dirvacuum==3)then
2979 1312 : rsuper(:,3)=rprimd(:,3)
2980 328 : shiftk_trial(3,1)=0.0_dp
2981 : end if
2982 :
2983 : ! The supercell and the corresponding shift have been generated !
2984 : ! Convert cartesian coordinates into kptrlatt_trial
2985 1368 : do ii=1,3
2986 : kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+&
2987 : & gprimd(2,:)*rsuper(2,ii)+&
2988 4446 : & gprimd(3,:)*rsuper(3,ii) )
2989 : end do
2990 :
2991 : ! End of 2-dimensional system
2992 : end if
2993 :
2994 : ! 3-dimensional system
2995 8216 : if(sum(vacuum(:))==0)then
2996 : ! Treat hexagonal holohedries separately
2997 1706 : if(iholohedry==6)then
2998 108 : length1=length_axis1*mult1
2999 108 : length3=length_axis3*mult3
3000 : ! write(std_out,*)' testkgrid: (hex) lengths=',length1,length2
3001 108 : if(abs(length1-length3)<tol8)then
3002 20 : mult1=mult1+1
3003 20 : mult3=mult3+1
3004 88 : else if(length1>length3)then
3005 44 : mult3=mult3+1
3006 44 : else if(length3>length1)then
3007 44 : mult1=mult1+1
3008 : end if
3009 108 : nset=4
3010 108 : if(iset==1)then
3011 100 : rsuper(:,1)=axes(:,1)*mult1
3012 100 : rsuper(:,2)=axes(:,2)*mult1
3013 100 : rsuper(:,3)=axes(:,3)*mult3
3014 100 : shiftk_trial(:,1)=0.0_dp
3015 25 : shiftk_trial(3,1)=0.5_dp
3016 83 : else if(iset==2)then
3017 116 : rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1
3018 116 : rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
3019 116 : rsuper(:,3)=axes(:,3)*mult3
3020 29 : shiftk_trial(1,1)=1.0_dp/3.0_dp
3021 29 : shiftk_trial(2,1)=1.0_dp/3.0_dp
3022 29 : shiftk_trial(3,1)=0.5_dp
3023 54 : else if(iset==3)then
3024 116 : rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1
3025 116 : rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
3026 116 : rsuper(:,3)=axes(:,3)*mult3
3027 116 : shiftk_trial(:,1)=0.0_dp
3028 29 : shiftk_trial(3,1)=0.5_dp
3029 25 : else if(iset==4)then
3030 100 : rsuper(:,1)=axes(:,1)*mult1
3031 100 : rsuper(:,2)=axes(:,2)*mult1
3032 100 : rsuper(:,3)=axes(:,3)*mult3
3033 100 : shiftk_trial(:,1)=0.5_dp
3034 : end if
3035 :
3036 : else
3037 : ! Now treat all other holohedries
3038 1598 : length1=length_axis1*mult1
3039 1598 : length2=length_axis2*mult2
3040 1598 : length3=length_axis3*mult3
3041 : ! write(std_out,*)' testkgrid: length=',length1,length2,length3
3042 1598 : if(length2>length1+tol8 .and. length3>length1+tol8)then
3043 100 : mult1=mult1+1
3044 1498 : else if(length1>length2+tol8 .and. length3>length2+tol8)then
3045 154 : mult2=mult2+1
3046 1344 : else if(length1>length3+tol8 .and. length2>length3+tol8)then
3047 352 : mult3=mult3+1
3048 : else if(abs(length2-length3)<tol8 .and. &
3049 992 : & abs(length1-length3)<tol8 .and. &
3050 : & abs(length1-length2)<tol8 )then
3051 807 : mult1=mult1+1 ; mult2=mult2+1 ; mult3=mult3+1
3052 185 : else if(abs(length1-length2)<tol8)then
3053 151 : mult1=mult1+1 ; mult2=mult2+1
3054 34 : else if(abs(length1-length3)<tol8)then
3055 0 : mult1=mult1+1 ; mult3=mult3+1
3056 34 : else if(abs(length2-length3)<tol8)then
3057 34 : mult2=mult2+1 ; mult3=mult3+1
3058 : end if
3059 1598 : nset=6
3060 1598 : if(center==-1 .or. center==-3)nset=8
3061 1598 : if(iset==1 .or. iset==2)then
3062 : ! Simple lattice of k points
3063 2080 : rsuper(:,1)=axes(:,1)*mult1
3064 2080 : rsuper(:,2)=axes(:,2)*mult2
3065 2080 : rsuper(:,3)=axes(:,3)*mult3
3066 2080 : shiftk_trial(:,1)=0.5_dp
3067 1078 : else if(iset==3 .or. iset==4)then
3068 : ! FCC lattice of k points = BCC lattice in real space
3069 1848 : rsuper(:,1)=-axes(:,1)*mult1+axes(:,2)*mult2+axes(:,3)*mult3
3070 1848 : rsuper(:,2)= axes(:,1)*mult1-axes(:,2)*mult2+axes(:,3)*mult3
3071 1848 : rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2-axes(:,3)*mult3
3072 1848 : shiftk_trial(:,1)=0.5_dp
3073 616 : else if(iset==5 .or. iset==6)then
3074 : ! BCC lattice of k points = FCC lattice in real space
3075 2028 : rsuper(:,1)= axes(:,2)*mult2+axes(:,3)*mult3
3076 2028 : rsuper(:,2)= axes(:,1)*mult1 +axes(:,3)*mult3
3077 2028 : rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2
3078 : ! The BCC lattice has no empty site with full symmetry
3079 2028 : shiftk_trial(:,1)=0.0_dp
3080 109 : else if(iset==7 .or. iset==8)then
3081 : ! iset==7 and 8 are allowed only for centered lattice
3082 109 : mult1h=mult1-0.5_dp
3083 109 : mult2h=mult2-0.5_dp
3084 109 : mult3h=mult3-0.5_dp
3085 109 : if(center==-1)then
3086 : ! FCC lattice of k points = BCC lattice in real space
3087 352 : rsuper(:,1)=-axes(:,1)*mult1h+axes(:,2)*mult2h+axes(:,3)*mult3h
3088 352 : rsuper(:,2)= axes(:,1)*mult1h-axes(:,2)*mult2h+axes(:,3)*mult3h
3089 352 : rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h-axes(:,3)*mult3h
3090 352 : shiftk_trial(:,1)=0.5_dp
3091 21 : else if(center==-3)then
3092 : ! BCC lattice of k points = FCC lattice in real space
3093 84 : rsuper(:,1)= axes(:,2)*mult2h+axes(:,3)*mult3h
3094 84 : rsuper(:,2)= axes(:,1)*mult1h +axes(:,3)*mult3h
3095 84 : rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h
3096 : ! The BCC lattice has no empty site with full symmetry
3097 84 : shiftk_trial(:,1)=0.0_dp
3098 : end if
3099 : end if
3100 : ! This was the easiest way to code all even mult1, mult2, mult3 triplets:
3101 : ! make separate series for this possibility.
3102 1598 : if(2*(iset/2)==iset)then
3103 2864 : rsuper(:,1)=2.0_dp*rsuper(:,1)
3104 2864 : rsuper(:,2)=2.0_dp*rsuper(:,2)
3105 2864 : rsuper(:,3)=2.0_dp*rsuper(:,3)
3106 : end if
3107 : end if
3108 :
3109 : ! write(std_out,*)' testkgrid: gprimd=',gprimd(:,:)
3110 : ! write(std_out,*)' testkgrid: rsuper=',rsuper(:,:)
3111 : ! write(std_out,*)' testkgrid: iset =',iset
3112 :
3113 : ! The supercell and the corresponding shift have been generated!
3114 : ! Convert cartesian coordinates into kptrlatt_trial
3115 6824 : do ii=1,3
3116 : kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+&
3117 : & gprimd(2,:)*rsuper(2,ii)+&
3118 22178 : & gprimd(3,:)*rsuper(3,ii) )
3119 : end do
3120 :
3121 : ! End of 3-dimensional system
3122 : end if
3123 :
3124 : ! write(std_out,*)' testkgrid: before getkgrid'
3125 : ! write(std_out,*)' testkgrid: rprimd=',rprimd(:,:)
3126 : ! write(std_out,*)' testkgrid: kptrlatt_trial=',kptrlatt_trial(:,:)
3127 :
3128 : call getkgrid(0,0,iscf,kpt,&
3129 : & kptopt,kptrlatt_trial,kptrlen_trial,&
3130 : & msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,&
3131 2054 : & shiftk_trial,symafm,symrel,vacuum,wtk)
3132 :
3133 : ! write(std_out,*)' testkgrid: after getkgrid'
3134 :
3135 : ! In case one does not need the full list of grids, will take a shortcut, and go to one of the last grids of the series,
3136 : ! that generates a kptrlen_trial that is just below kptrlen.
3137 2054 : if(prtkpt==0 .and. init_mult==1 .and. kptrlen_trial<(half-tol8)*kptrlen )then
3138 274 : iscale=int((one-tol8)*kptrlen/kptrlen_trial)
3139 274 : mult1=mult1*iscale
3140 274 : mult2=mult2*iscale
3141 274 : mult3=mult3*iscale
3142 274 : init_mult=0
3143 : ! write(std_out,*)' testkgrid: iscale=',iscale
3144 3562 : kptrlatt_trial(:,:)=kptrlatt_trial(:,:)*iscale
3145 : call getkgrid(0,0,iscf,kpt,&
3146 : & kptopt,kptrlatt_trial,kptrlen_trial,&
3147 : & msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,&
3148 274 : & shiftk_trial,symafm,symrel,vacuum,wtk)
3149 : end if
3150 :
3151 : if( (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_current==0) .or. &
3152 2054 : & (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_trial<nkpt_current) .or. &
3153 : & (nkpt_trial==nkpt_current .and. kptrlen_trial>kptrlen_current*(1.0_dp+tol8)))then
3154 :
3155 181 : kptrlatt_current(:,:)=kptrlatt_trial(:,:)
3156 181 : nkpt_current=nkpt_trial
3157 181 : shiftk_current(:,:)=shiftk_trial(:,:)
3158 181 : kptrlen_current=kptrlen_trial
3159 181 : igrid_current=igrid
3160 : end if
3161 :
3162 2054 : if(prtkpt/=0)then
3163 : write(msg,'(i5,a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )&
3164 788 : & igrid,' ',kptrlatt_trial(:,1),' ',shiftk_trial(1,1),&
3165 788 : & ' ',kptrlen_trial,nkpt_trial,iset,ch10,&
3166 788 : & ' ',kptrlatt_trial(:,2),' ',shiftk_trial(2,1),ch10,&
3167 1576 : & ' ',kptrlatt_trial(:,3),' ',shiftk_trial(3,1),ch10
3168 788 : call wrtout(std_out,msg)
3169 788 : call wrtout(iout,msg)
3170 :
3171 : ! Keep track of this grid, if it is worth
3172 788 : if(kptrlen_trial > kptrlen_list(nkpt_trial)*(1.0_dp+tol8))then
3173 512 : grid_list(nkpt_trial)=igrid
3174 512 : kptrlen_list(nkpt_trial)=kptrlen_trial
3175 : end if
3176 : end if
3177 :
3178 : ! Treat 1-D case
3179 8216 : if( sum(vacuum(:))==2 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )exit
3180 :
3181 : ! Treat 2-D case or 3-D case
3182 8295 : if( sum(vacuum(:))<=1 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )then
3183 : ! The present set of sets of k points is finished:
3184 : ! either it was the last, or one has to go to the next one
3185 564 : if(iset==nset)exit
3186 473 : iset=iset+1
3187 473 : mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1
3188 : end if
3189 :
3190 : end do ! igrid=1,1000
3191 :
3192 95 : ABI_FREE(kpt)
3193 95 : ABI_FREE(wtk)
3194 :
3195 95 : kptrlatt(:,:)=kptrlatt_current(:,:)
3196 95 : shiftk(:,:)=shiftk_current(:,:)
3197 95 : kptrlen=kptrlen_current
3198 :
3199 : end if ! test on the number of dimensions
3200 :
3201 197 : if(prtkpt/=0)then
3202 :
3203 : ! sqrt(1/2) comes from the FCC packing, the best one
3204 13 : factor=sqrt(0.5_dp)/ucvol/dble(nsym)
3205 13 : ndims=3
3206 52 : if(sum(vacuum(:))/=0)then
3207 24 : if(sum(vacuum(:))==1)then
3208 : ! sqrt(3/4) comes from the hex packing, the best one
3209 : ! one multiplies by 2 because nsym is likely twice the number
3210 : ! of symmetries that can be effectively used in 2D
3211 6 : ndims=2 ; factor=sqrt(0.75_dp)/surface/dble(nsym)*2
3212 6 : write(msg,'(2a)' )ch10,' Note that the system is bi-dimensional.'
3213 0 : else if(sum(vacuum(:))==2)then
3214 0 : ndims=1 ; factor=1/ucvol
3215 0 : write(msg,'(2a)' )ch10,' Note that the system is uni-dimensional.'
3216 0 : else if(sum(vacuum(:))==3)then
3217 0 : ndims=0
3218 0 : write(msg,'(2a)' )ch10,' Note that the system is zero-dimensional.'
3219 : end if
3220 6 : call wrtout(std_out,msg)
3221 6 : call wrtout(iout,msg)
3222 : end if
3223 :
3224 : ! The asymptotic value of the merit factor is determined
3225 : ! by the set of symmetries: in 3D, if it includes the
3226 : ! inversion symmetry, the limit will be 1, if not, it
3227 : ! will be two. In 2D, if it includes the inversion symmetry
3228 : ! and an operation that maps z on -z, it will tend to one,
3229 : ! while if only one of these operations is present,
3230 : ! it will tend to two, and if none is present, it will tend to four.
3231 13 : write(msg,'(11a)' )ch10,&
3232 13 : & ' List of best grids, ordered by nkpt.',ch10,&
3233 13 : & ' (stop at a value of kptrlen 20% larger than the target value).',ch10,&
3234 13 : & ' (the merit factor will tend to one or two in 3 dimensions)',ch10,&
3235 13 : & ' (and to one, two or four in 2 dimensions)',ch10,ch10,&
3236 26 : & ' nkpt kptrlen grid# merit_factor'
3237 13 : call wrtout(std_out,msg)
3238 13 : call wrtout(iout,msg)
3239 :
3240 13 : kptrlen_max=0.0_dp
3241 347 : do ii=1,mkpt_list
3242 347 : if(kptrlen_list(ii)>kptrlen_max*(1.0_dp+tol8))then
3243 126 : kptrlen_max=kptrlen_list(ii)
3244 126 : merit_factor=kptrlen_max**ndims/dble(ii)*factor
3245 126 : write(msg, '(i6,es14.4,i6,f12.4)' )ii,kptrlen_max,grid_list(ii),merit_factor
3246 126 : call wrtout(std_out,msg)
3247 126 : call wrtout(iout,msg)
3248 : end if
3249 347 : if(kptrlen_max>1.2_dp*(1.0_dp-tol8)*kptrlen_target)exit
3250 : end do
3251 :
3252 13 : write(msg,'(a,a,es14.4,a,a,i6,a,a,a,es14.4,a,i6)' )ch10,&
3253 13 : & ' For target kptrlen=',kptrlen_target,',',&
3254 13 : & ' the selected grid is number',igrid_current,',',ch10,&
3255 26 : & ' giving kptrlen=',kptrlen_current,' with nkpt=',nkpt_current
3256 13 : call wrtout(std_out,msg)
3257 13 : call wrtout(iout,msg)
3258 :
3259 13 : write(msg,'(a,a,a,a)' )ch10,&
3260 13 : & ' testkgrid : stop after analysis of a series of k-grids.',ch10,&
3261 26 : & ' For usual production runs, set prtkpt back to 0 (the default).'
3262 13 : call wrtout(std_out,msg, do_flush=.True.)
3263 13 : call wrtout(iout,msg, do_flush=.True.)
3264 :
3265 13 : call abi_abort('PERS', exit_status=0, print_config=.false.)
3266 : end if
3267 :
3268 184 : end subroutine testkgrid
3269 : !!***
3270 :
3271 : !!****f* m_kpts/mknormpath
3272 : !! NAME
3273 : !! mknormpath
3274 : !!
3275 : !! FUNCTION
3276 : !! Please do not use this routine, use make_normpath instead.
3277 : !! mknormpath should be removed
3278 : !!
3279 : !! This simple routine generates a normalized path that can be used to plot a band
3280 : !! structures in an easy way. For normalized path we mean a path where the number
3281 : !! of division on each segment is proportional to the length of the segment itself.
3282 : !! To generate the above mentioned path, the subroutine must be called twice.
3283 : !! The first call reports the total number of divisions in the normalized path, dimension
3284 : !! that is required to correctly allocate the array.
3285 : !! The second call calculates the reduced coordinates of the circuit.
3286 : !!
3287 : !! INPUTS
3288 : !! nbounds=number of points defining the path
3289 : !! ndiv_small=number of points to be used to sample the smallest
3290 : !! segment defined by bounds(:,1:nbounds)
3291 : !! bounds(3,nbounds)=points defining the path
3292 : !! gmet(3,3)=metric
3293 : !!
3294 : !! OUTPUT
3295 : !! ndiv(nbounds-1)= number of divisions for each segment
3296 : !! npt_tot=total number of points sampled along the circuit
3297 : !! path(3,npt_tot)= normalized path in reciprocal space
3298 : !!
3299 : !! TODO
3300 : !! Do not use this routine, it is obsolete and should be replaced by make_path in m_bz_mesh.
3301 : !!
3302 : !! SOURCE
3303 :
3304 70 : subroutine mknormpath(nbounds,bounds,gmet,ndiv_small,ndiv,npt_tot,path)
3305 :
3306 : !Arguments ------------------------------------
3307 : !scalars
3308 : integer,intent(in) :: nbounds,ndiv_small
3309 : integer,intent(inout) :: npt_tot
3310 : !arrays
3311 : integer,intent(inout) :: ndiv(nbounds-1)
3312 : real(dp),intent(in) :: bounds(3,nbounds),gmet(3,3)
3313 : real(dp),intent(out),optional :: path(3,npt_tot)
3314 :
3315 : !Local variables-------------------------------
3316 : !scalars
3317 : integer :: idx,ii,jp
3318 : real(dp) :: fct
3319 : character(len=500) :: msg
3320 : !arrays
3321 140 : real(dp) :: dd(3),lng(nbounds-1)
3322 : ! *************************************************************************
3323 :
3324 70 : if (ndiv_small<=0) then
3325 : write(msg,'(3a,i0)')&
3326 0 : 'The argument ndiv_small should be a positive number,',ch10,&
3327 0 : 'however, ndiv_small=',ndiv_small
3328 0 : ABI_ERROR(msg)
3329 : end if
3330 :
3331 586 : do ii=1,nbounds-1
3332 2064 : dd(:)=bounds(:,ii+1)-bounds(:,ii)
3333 : lng(ii)= sqrt( dd(1)*gmet(1,1)*dd(1)+ &
3334 : & dd(2)*gmet(2,2)*dd(2)+ &
3335 : & dd(3)*gmet(3,3)*dd(3)+ &
3336 : & 2.0d0*(dd(1)*gmet(1,2)*dd(2)+ &
3337 : & dd(1)*gmet(1,3)*dd(3)+ &
3338 : & dd(2)*gmet(2,3)*dd(3)) &
3339 586 : & )
3340 : end do
3341 70 : write(std_out,*)lng
3342 656 : fct=minval(lng)
3343 :
3344 : ! Avoid division by zero if k(:,i+1)=k(:,i)
3345 70 : if (abs(fct)<tol6) then
3346 : write(msg,'(3a)')&
3347 0 : 'found two consecutive points in the path which are equal',ch10,&
3348 0 : 'This is not allowed, please modify the path in your input file'
3349 0 : ABI_ERROR(msg)
3350 : end if
3351 :
3352 70 : fct=fct/ndiv_small
3353 586 : ndiv(:)=nint(lng(:)/fct)
3354 : ! The 1 stand for the first point
3355 586 : npt_tot=sum(ndiv)+1
3356 :
3357 70 : if (.not.present(path)) then
3358 70 : write(msg,'(2a,i0)')ch10,' mknormpath : total number of points on the path: ',npt_tot
3359 70 : call wrtout(std_out,msg)
3360 70 : write(msg,'(2a)')ch10,' Number of divisions for each segment of the normalized path: '
3361 70 : call wrtout(std_out,msg)
3362 586 : do ii=1,nbounds-1
3363 : write(msg,'(2(3f8.5,a),i5,a)')&
3364 516 : bounds(:,ii),' ==> ',bounds(:,ii+1),' ( ndiv: ',ndiv(ii),' )'
3365 586 : call wrtout(std_out,msg)
3366 : end do
3367 70 : write(msg,'(a)')ch10
3368 70 : call wrtout(std_out,msg)
3369 : else
3370 0 : write(msg,'(2a)')ch10,' Normalized Path: '
3371 0 : call wrtout(std_out,msg)
3372 0 : idx=1
3373 0 : do ii=1,nbounds-1
3374 0 : do jp=1,ndiv(ii)
3375 0 : path(:,idx)=bounds(:,ii)+(jp-1)*(path(:,ii+1)-path(:,ii))/ndiv(ii)
3376 0 : write(msg,'(i4,4x,3(f8.5,1x))')idx,path(:,idx)
3377 0 : call wrtout(std_out,msg)
3378 0 : idx=idx+1
3379 : end do
3380 : end do
3381 : end if
3382 :
3383 70 : end subroutine mknormpath
3384 : !!***
3385 :
3386 : !!****f* m_kpts/bzlint_init
3387 : !! NAME
3388 : !! bzlint_init
3389 : !!
3390 : !! FUNCTION
3391 : !!
3392 : !! INPUTS
3393 : !!
3394 : !! OUTPUT
3395 : !!
3396 : !! SOURCE
3397 :
3398 0 : subroutine bzlint_init(self, ngkpt, ndat, nkpt, kpts, values)
3399 :
3400 : !Arguments ------------------------------------
3401 : class(bzlint_t),intent(inout) :: self
3402 : integer,intent(in) :: ngkpt(3), ndat, nkpt
3403 : real(dp),intent(in) :: kpts(3,nkpt), values(ndat, nkpt)
3404 :
3405 : !Local variables-------------------------------
3406 : integer :: ik, ix, iy, iz, inds(3)
3407 : real(dp) :: kpt_wrap(3), shift(3)
3408 : ! *********************************************************************
3409 :
3410 0 : self%ngkpt = ngkpt; self%ndat = ndat
3411 : ! The mesh is closed i.e. periodic images are included.
3412 0 : self%nx = ngkpt(1); self%ny = ngkpt(2); self%nz = ngkpt(3)
3413 0 : ABI_CALLOC(self%vals_grid, (self%nx, self%ny, self%nz, ndat))
3414 :
3415 : ! Insert values in the grid.
3416 0 : do ik=1,nkpt
3417 0 : call wrap2_zero_one(kpts(:,ik), kpt_wrap, shift)
3418 0 : inds = nint(kpt_wrap * self%ngkpt)
3419 : ! here we need to shift the indices by 1 (since Fortran people like to count starting from 1 (: )
3420 0 : ix = inds(1) + 1; iy = inds(2) + 1; iz = inds(3) + 1
3421 0 : self%vals_grid(ix,iy,iz,:) = values(:, ik)
3422 : end do
3423 :
3424 0 : end subroutine bzlint_init
3425 : !!***
3426 :
3427 : !!****f* m_kpts/bzlint_interp
3428 : !! NAME
3429 : !! bzlint_interp
3430 : !!
3431 : !! FUNCTION
3432 : !! Interpolate values at kpt
3433 : !!
3434 : !! SOURCE
3435 :
3436 0 : subroutine bzlint_interp(self, kpt, results)
3437 :
3438 : !Arguments ------------------------------------
3439 : class(bzlint_t),intent(in) :: self
3440 : real(dp),intent(in) :: kpt(3)
3441 : real(dp),intent(out) :: results(self%ndat)
3442 :
3443 : !Local variables-------------------------------
3444 : integer :: idat
3445 : real(dp) :: kpt_wrap(3), shift(3)
3446 : ! *********************************************************************
3447 :
3448 0 : call wrap2_zero_one(kpt, kpt_wrap, shift)
3449 0 : do idat=1,self%ndat
3450 0 : results(idat) = interpol3d_0d(kpt_wrap, self%nx, self%ny, self%nz, self%vals_grid(:,:,:,idat))
3451 : end do
3452 :
3453 0 : end subroutine bzlint_interp
3454 : !!***
3455 :
3456 : !!****f* m_kpts/bzlint_free
3457 : !! NAME
3458 : !! bzlint_free
3459 : !!
3460 : !! FUNCTION
3461 : !! Free dynamic memory
3462 : !!
3463 : !! SOURCE
3464 :
3465 0 : subroutine bzlint_free(self)
3466 :
3467 : !Arguments ------------------------------------
3468 : class(bzlint_t),intent(inout) :: self
3469 : ! *********************************************************************
3470 :
3471 0 : ABI_SFREE(self%vals_grid)
3472 :
3473 0 : end subroutine bzlint_free
3474 : !!***
3475 :
3476 : !!****f* m_kpts/kptrlatt_from_ngkpt
3477 : !! NAME
3478 : !! kptrlatt_from_ngkpt
3479 : !!
3480 : !! FUNCTION
3481 : !! Insert ngkpt in kptrlatt 3x3 array
3482 : !!
3483 : !! SOURCE
3484 :
3485 4 : pure subroutine kptrlatt_from_ngkpt(ngkpt, kptrlatt)
3486 :
3487 : !Arguments ------------------------------------
3488 : integer,intent(in) :: ngkpt(3)
3489 : integer,intent(out) :: kptrlatt(3,3)
3490 :
3491 : !Local variables-------------------------------
3492 : integer :: ii
3493 : !************************************************************************
3494 :
3495 4 : kptrlatt = 0
3496 16 : do ii=1,3
3497 16 : kptrlatt(ii,ii) = ngkpt(ii)
3498 : end do
3499 :
3500 4 : end subroutine kptrlatt_from_ngkpt
3501 : !!***
3502 :
3503 0 : end module m_kpts
3504 : !!***
|
