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