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