Line data Source code
1 : !!****m* ABINIT/m_bz_mesh
2 : !! NAME
3 : !! m_bz_mesh
4 : !!
5 : !! FUNCTION
6 : !! This module provides the definition of the kmesh_t structure gathering information
7 : !! on the sampling of the Brillouin zone. It also contains useful tools to operate on k-points.
8 : !! and the definition of the littlegroup_t data type. The littlegroup_t structure is used
9 : !! to store tables and useful info on the set of k-points belonging
10 : !! to the irreducible wedge defined by the symmetry properties
11 : !! of the point group that preserve the external q-point.
12 : !!
13 : !! COPYRIGHT
14 : !! Copyright (C) 2008-2026 ABINIT group (MG, GMR, VO, LR, RWG, MT)
15 : !! This file is distributed under the terms of the
16 : !! GNU General Public License, see ~abinit/COPYING
17 : !! or http://www.gnu.org/copyleft/gpl.txt .
18 : !!
19 : !! NOTES
20 : !! One has to use a fixed ordering of the loops over nsym and time-reversal
21 : !! when the full zone is reconstructed by symmetry starting from the IBZ.
22 : !! This is especially important in systems with both time-reversal and
23 : !! spatial inversion as the effect of the two operation in reciprocal
24 : !! space is very similar the only difference being the possibly non-zero
25 : !! fractional translation associated to the spatial inversion.
26 : !! In the present implementation, the spatial inversion has the precedence
27 : !! wrt time-reversal (i.e., do itim; do isym).
28 : !! Note that this particular ordering should be used in any routine used to
29 : !! symmetrize k-dependent quantities in the full BZ zone to avoid possible errors.
30 : !!
31 : !! Important WARNING:
32 : !!
33 : !! This module is deprecated and should be used only in the GW/BSE part.
34 : !! Some of the routines will be gradually moved to m_kpts
35 : !!
36 : !! SOURCE
37 :
38 : #if defined HAVE_CONFIG_H
39 : #include "config.h"
40 : #endif
41 :
42 : #include "abi_common.h"
43 :
44 : module m_bz_mesh
45 :
46 : use defs_basis
47 : use m_errors
48 : use m_abicore
49 : use m_sort
50 : use m_xmpi
51 :
52 : use m_fstrings, only : ltoa, itoa, sjoin, ktoa
53 : use m_numeric_tools, only : is_zero, isinteger, imin_loc, imax_loc, bisect, wrap2_pmhalf
54 : use m_symtk, only : sg_multable, littlegroup_q
55 : use m_geometry, only : normv
56 : use m_crystal, only : crystal_t
57 : use m_kpts, only : getkgrid
58 : use m_symkpt, only : symkpt
59 :
60 : implicit none
61 :
62 : private
63 :
64 : real(dp),parameter :: TOL_KDIFF = 0.0001_dp
65 : ! Tolerance below which two points are considered equal within a RL vector:
66 : ! for each reduced direction the absolute difference between the coordinates must be less that TOL_KDIFF
67 :
68 : integer,parameter :: NONE_KPTRLATT(3,3) = RESHAPE((/0,0,0,0,0,0,0,0,0/),(/3,3/))
69 : !!***
70 :
71 : !!****t* m_bz_mesh/kmesh_t
72 : !! NAME
73 : !! kmesh_t
74 : !!
75 : !! FUNCTION
76 : !! The kmesh_t structured datatype contains different information on the grid used to sample the BZ :
77 : !! the k-points in the full Brillouin zone BZ, the irreducible wedge IBZ as well as tables describing
78 : !! the symmetry relationship between the points.
79 : !!
80 : !! SOURCE
81 :
82 : type,public :: kmesh_t
83 :
84 : !scalars
85 : integer :: nshift = 0
86 :
87 : integer :: nbz = 0
88 : ! Number of points in the BZ.
89 :
90 : integer :: nibz = 0
91 : ! Number of points in the IBZ.
92 :
93 : integer :: nsym
94 : ! Number of symmetry operations.
95 :
96 : integer :: kptopt
97 : ! Option for the generation of k points (see input variable description)
98 : !
99 : ! 1 if both time-reversal and point group symmetries are used.
100 : ! 2 if only time-reversal symmetry is used.
101 : ! 3 do not take into account any symmetry (except the identity).
102 : ! 4 if time-reversal is not used (spin-orbit coupling).
103 : ! < 0 number of segments used to construct the k-path for NSCF calculation.
104 :
105 : integer :: timrev
106 : ! 2 if time reversal symmetry can be used, 1 otherwise.
107 :
108 : !arrays
109 : integer :: kptrlatt(3,3) = NONE_KPTRLATT
110 : ! Coordinates of three vectors in real space, expressed in reduced coordinates.
111 : ! They define a super-lattice in real space. The k point lattice is the reciprocal of
112 : ! this super-lattice, eventually shifted by shift.
113 : ! Not available if the structure is initialized from the points in the IBZ.
114 :
115 : integer,allocatable :: rottb(:,:,:)
116 : ! (nbz, timrev, nsym),
117 : ! Index of (IS)k in the BZ array where S is a sym operation in reciprocal space,
118 : ! I is the identity or the inversion operator (1,2 resp)
119 :
120 : integer,allocatable :: rottbm1(:,:,:)
121 : ! (nbz, timrev, nsym)
122 : ! Index of IS^{-1} k in the BZ array.
123 :
124 : integer,allocatable :: tab(:)
125 : ! (nbz)
126 : ! For each point in the BZ, it gives the index of the symmetric irreducible point in the ibz.
127 :
128 : integer,allocatable :: tabi(:)
129 : ! (nbz)
130 : ! For each point in the BZ, tabi tells whether time-reversal has to be
131 : ! used to obtain k_BZ starting from the corresponding point in the IBZ (1=>no, -1=>yes)
132 :
133 : integer,allocatable :: tabo(:)
134 : ! (nbz)
135 : ! For each point in the BZ, it gives the index in the array symrec of the
136 : ! symmetry operation in reciprocal space which rotates k_IBZ onto \pm k_BZ (depending on tabi)
137 :
138 : integer,allocatable :: umklp(:,:)
139 : ! (3,nbz)
140 : ! The Umklapp G0-vector such as kbz + G0 = (IS) k_ibz, where kbz is in the first BZ.
141 :
142 : real(dp) :: gmet(3,3)
143 : ! Reciprocal space metric ($\textrm{bohr}^{-2}$).
144 :
145 : real(dp) :: gprimd(3,3)
146 : ! Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
147 :
148 : real(dp),allocatable :: bz(:,:)
149 : ! (3,nbz)
150 : ! Points in the BZ in reduced coordinates.
151 : ! TODO should packed in shells.
152 :
153 : real(dp),allocatable :: ibz(:,:)
154 : ! (3, nibz)
155 : ! Points in the IBZ in reduced coordinates.
156 :
157 : real(dp),allocatable :: shift(:,:)
158 : ! (3, nshift)
159 : ! Shift for k-points, not available is nshift=0. Usually nshift=1
160 :
161 : real(dp),allocatable :: wt(:)
162 : ! (nibz)
163 : ! Weights for each point in the IBZ.
164 :
165 : complex(dp),allocatable :: tabp(:)
166 : ! (nkbz)
167 : ! For each point in the BZ, this table gives the phase factors associated
168 : ! to non-symmorphic operations, i.e., e^{-i2\pi k_IBZ.R{^-1}t}=e^{-i2\pi k_BZ cdot t}
169 : ! where \transpose R{-1}=S and (S k_IBZ)=\pm k_BZ (depending on tabi)
170 :
171 : contains
172 :
173 : ! Methods
174 : procedure :: init => kmesh_init ! Main creation method.
175 : procedure :: free => kmesh_free ! Free memory
176 : procedure :: print => kmesh_print ! Printout of basic info on the object.
177 : procedure :: get_bz_item => get_bz_item ! Get point in the BZ and other useful quantities.
178 : procedure :: get_ibz_item => get_IBZ_item ! Get point in the IBZ and other useful quantities.
179 : procedure :: get_bz_diff => get_BZ_diff ! Get the difference k1-k2 in the BZ (if any).
180 : procedure :: has_bz_item => has_BZ_item ! Check if a point belongs to the BZ mesh.
181 : procedure :: has_ibz_item => has_IBZ_item ! Check if a point is in the IBZ
182 : procedure :: isirred => bz_mesh_isirred ! TRUE if ik_bz is in the IBZ (a non-zero umklapp is not allowed)
183 : !prodedure :: pack_in_stars => bz_mesh_pack_in_stars
184 :
185 : procedure :: find_qmesh => bz_mesh_find_qmesh ! Find the Q-mesh defined as the set of all possible k1-k2 differences.
186 :
187 : end type kmesh_t
188 :
189 : public :: make_mesh ! Initialize the mesh starting from kptrlatt and shiftk.
190 : public :: isamek ! Check whether two points are equal within an umklapp G0.
191 : public :: isequalk ! Check whether two points are equal within an umklapp G0 (does not report G0)
192 : public :: findqg0 ! Identify q + G0 = k1-k2.
193 : public :: findq ! Helper routine returning the list of q-points.
194 : public :: findnq ! Helper routine returning the number of q-points.
195 : public :: identk ! Find the BZ starting from the irreducible k-points.
196 : public :: get_ng0sh ! Calculate the smallest box in RSpace able to treat all possible umklapp processes.
197 : public :: box_len ! Return the length of the vector connecting the origin with one the faces of the unit cell.
198 : !!***
199 :
200 : !----------------------------------------------------------------------
201 :
202 : !!****t* m_bz_mesh/kpath_t
203 : !! NAME
204 : !! path_t
205 : !!
206 : !! FUNCTION
207 : !! A (normalized) path in reciprocal space
208 : !!
209 : !! SOURCE
210 :
211 : type,public :: kpath_t
212 :
213 : integer :: nbounds = 0
214 : ! Number of extrema defining the path.
215 :
216 : integer :: ndivsm = 0
217 : ! Number of divisions used to sample the smallest segment.
218 :
219 : integer :: npts = 0
220 : ! Total number of points in the path.
221 :
222 : real(dp) :: gprimd(3,3)
223 : ! Reciprocal lattice vectors.
224 :
225 : real(dp) :: gmet(3,3)
226 : ! Metric matrix in G-space.
227 :
228 : integer,allocatable :: ndivs(:)
229 : ! ndivs(nbounds-1)
230 : ! Number of divisions for each segment.
231 :
232 : integer,allocatable :: bounds2kpt(:)
233 : ! (nbounds)
234 : ! bounds2kpt(i): Index of the i-th extrema in the pts(:) array.
235 :
236 : real(dp),allocatable :: bounds(:,:)
237 : ! (3,nbounds)
238 : ! The points defining the path in reduced coordinates.
239 :
240 : real(dp),allocatable :: points(:,:)
241 : ! (3,npts)
242 : ! The points of the path in reduced coordinates.
243 :
244 : real(dp),allocatable :: dl(:)
245 : ! (npts)
246 : ! dl(i) = Distance between the (i-1)-th and the i-th k-point. dl(1) = zero
247 :
248 : contains
249 :
250 : procedure :: init => kpath_init
251 : ! Construct a new path
252 :
253 : procedure :: free => kpath_free
254 : ! Free memory
255 :
256 : procedure :: get_versors => kpath_get_versors
257 : ! Return all the versors emanating from the Gamma point.
258 :
259 : procedure :: print => kpath_print
260 : ! Print the path.
261 :
262 : end type kpath_t
263 :
264 :
265 : public :: make_path ! Construct a normalized path. TODO: Remove it as it's deprecated
266 : !!***
267 :
268 : !----------------------------------------------------------------------
269 :
270 : !!****t* m_bz_mesh/littlegroup_t
271 : !! NAME
272 : !! littlegroup_t
273 : !!
274 : !! FUNCTION
275 : !! For the GW part of ABINIT. The littlegroup_t structured datatype gather information on
276 : !! the little group associated to an external vector q. The little group associated to q
277 : !! is defined as the subset of the space group that preserves q, modulo a G0 vector
278 : !! (also called umklapp vector). Namely
279 : !!
280 : !! Sq = q +G0, where S is an operation in reciprocal space.
281 : !!
282 : !! If time reversal symmetry holds true, it is possible to enlarge the little group by
283 : !! including the operations such as
284 : !! -Sq = q+ G0.
285 : !!
286 : !! The operations belonging to the little group define an irriducible wedge in the Brillouin zone
287 : !! that is, usually, larger than the irredubile zone defined by the space group.
288 : !! The two zone coincide when q=0
289 : !!
290 : !! TODO
291 : !! Rationalize most of the arrays, in particular the tables
292 : !! This structure should be rewritten almost from scratch, thus avoid using it
293 : !! for your developments.
294 : !!
295 : !! SOURCE
296 :
297 : type,public :: littlegroup_t
298 :
299 : integer :: npw ! No. of planewaves used to describe the wavefuntion, used to dimension igmG0
300 : integer :: nsym_sg ! No. of operations in the space group (*NOT* the little group)
301 : integer :: nsym_ltg ! No. of symmetry operations in the little group (time-reversal is included, if can be used)
302 : integer :: timrev ! 2 if time-reversal is considered, 1 otherwise
303 : integer :: nbz ! No. of kpoints in the full BZ
304 : integer :: nibz_ltg ! No. of points in the irreducible wedge defined by the little group
305 : !integer :: use_umklp ! 1 if umklapp processes are included
306 :
307 : real(dp) :: max_kin_gmG0
308 : ! Max kinetic energy of G-G0 in case of umklapp.
309 :
310 : integer,allocatable :: G0(:,:,:)
311 : ! (3,2,nsym_sg)
312 : ! Reduced coordinates of the umklapp G0 vector.
313 :
314 : integer,allocatable :: ibzq(:)
315 : ! ibzq(nbz)
316 : ! 1 if the point belongs to the IBZ_q defined by ext_pt, 0 otherwise.
317 :
318 : integer,allocatable :: bz2ibz(:)
319 : ! bz2ibz(nbz)
320 : ! Index of the point in the irreducible wedge defined by the little group, 0 otherwise.
321 :
322 : integer,allocatable :: ibz2bz(:)
323 : ! ibz2bz(nibz_ltg)
324 : ! The corresponding index in the BZ array
325 :
326 : integer,allocatable :: igmG0(:,:,:)
327 : ! iumklp(npw,2,nsym_sg)
328 : ! Index of G-G0 in the FFT array for each operations IS (I=\pm 1).
329 :
330 : integer,allocatable :: flag_umklp(:,:)
331 : ! flag_umklp(2,nsym_sg)
332 : ! 1 if the operation IS requires a non null G0 vector to preserve q, 0 otherwise.
333 :
334 : integer,allocatable :: preserve(:,:)
335 : ! preserve(2, nsym_sg)
336 : ! (1,S) is 1 if the operation S in rec space preserves the external q-point i.e Sq=q+G0
337 : ! (2,S) is 1 if -Sq=q+G0. G0 is a reciprocal lattice vector also called "umklapp vector".
338 :
339 : integer,allocatable :: tab(:)
340 : ! tab(nbz)
341 : ! For each point in BZ, the index of the irreducible point (kIBZ_q) in the irreducible
342 : ! wedge defined by the little group of q. kBZ= (IS) kIBZ where I is the inversion or the identity.
343 :
344 : integer,allocatable :: tabo(:)
345 : ! tabo(nbz)
346 : ! The index of the operation S in the little group that rotates kIBZ_q into \pm kBZ.
347 :
348 : integer,allocatable :: tabi(:)
349 : ! tabi(nbz)
350 : ! for each k-point in the BZ defines whether inversion has to be
351 : ! considered in the relation kBZ= IS kIBZ_q (1 => only S; -1 => -S).
352 :
353 : integer,allocatable :: wtksym(:,:,:)
354 : ! (2, nsym_sg, kbz)
355 : ! 1 if IS belongs to the little group, 0 otherwise TODO (should invert the first two dimensions)
356 :
357 : real(dp) :: ext_pt(3)
358 : ! The external point defining the little group.
359 :
360 : contains
361 :
362 : procedure :: init => littlegroup_init
363 : procedure :: print => littlegroup_print
364 : procedure :: free => littlegroup_free_0D
365 : !procedure :: find => littlegroup_find
366 :
367 : end type littlegroup_t
368 :
369 : public :: littlegroup_free
370 : !!***
371 :
372 : interface littlegroup_free
373 : module procedure littlegroup_free_0D
374 : module procedure littlegroup_free_1D
375 : end interface littlegroup_free
376 :
377 : CONTAINS !=============================================================================
378 : !!***
379 :
380 : !!****f* m_bz_mesh/kmesh_init
381 : !! NAME
382 : !! kmesh_init
383 : !!
384 : !! FUNCTION
385 : !! Initialize and construct a kmesh_t datatype
386 : !! gathering information on the mesh in the Brilloin zone.
387 : !!
388 : !! INPUTS
389 : !! nkibz=Number of irreducible k-points.
390 : !! kibz(3,nkibz)=Irreducible k-points in reduced coordinates.
391 : !! Cryst<crystal_t> = Info on unit cell and its symmetries
392 : !! %nsym=number of symmetry operations
393 : !! %symrec(3,3,nsym)=symmetry operations in reciprocal space
394 : !! %tnons(3,nsym)=fractional translations
395 : !! kptopt=option for the generation of k points (see input variable description)
396 : !! [wrap_1zone]=If .TRUE., the points are wrapped in in the first BZ. Defaults to .FALSE. to preserve GW implementation.
397 : !! [ref_bz(:,:)]= Reference set of points in the full Brillouin zone used to prune k-points.
398 : !!
399 : !! OUTPUT
400 : !! Kmesh<kmesh_t>=Datatype gathering information on the k point sampling.
401 : !!
402 : !! SOURCE
403 :
404 1612 : subroutine kmesh_init(Kmesh, cryst, nkibz, kibz, kptopt, &
405 806 : wrap_1zone, ref_bz, break_symmetry) ! Optional
406 :
407 : !Arguments ------------------------------------
408 : !scalars
409 : class(kmesh_t),intent(inout) :: Kmesh
410 : integer,intent(in) :: nkibz, kptopt
411 : logical,optional,intent(in) :: wrap_1zone,break_symmetry
412 : type(crystal_t),intent(in) :: Cryst
413 : !arrays
414 : real(dp),intent(in) :: kibz(3,nkibz)
415 : real(dp),optional,intent(in) :: ref_bz(:,:)
416 :
417 : !Local variables-------------------------------
418 : !scalars
419 : integer :: ik_bz,ik_ibz,isym,nkbz,nkbzX,nsym,timrev,itim
420 : real(dp) :: shift(3)
421 : logical :: ltest,do_wrap,do_hack
422 : !arrays
423 806 : integer,allocatable :: ktab(:),ktabi(:),ktabo(:)
424 : real(dp) :: rm1t(3),kbz_wrap(3)
425 806 : real(dp),allocatable :: kbz(:,:),wtk(:)
426 : ! *************************************************************************
427 :
428 : ! === Initial tests on input arguments ===
429 806 : ltest = (Cryst%timrev==1 .or. Cryst%timrev==2)
430 806 : ABI_CHECK(ltest, sjoin('Wrong value for timrev= ', itoa(Cryst%timrev)))
431 :
432 806 : if (all(kptopt/= [1, 3])) then
433 14 : ABI_WARNING(sjoin("Not allowed value for kptopt: ", itoa(kptopt)))
434 : end if
435 :
436 806 : Kmesh%kptopt = kptopt
437 806 : nsym = Cryst%nsym; timrev = Cryst%timrev
438 :
439 : ! Find BZ from IBZ and fill tables ===
440 806 : nkbzX = nkibz*nsym*timrev ! Maximum possible number
441 2418 : ABI_MALLOC(kbz, (3,nkbzX))
442 2418 : ABI_MALLOC(wtk, (nkibz))
443 2418 : ABI_MALLOC(ktab, (nkbzX))
444 1612 : ABI_MALLOC(ktabi, (nkbzX))
445 1612 : ABI_MALLOC(ktabo, (nkbzX))
446 :
447 806 : if (present(ref_bz)) then
448 22 : call identk(kibz,nkibz,nkbzX,nsym,timrev,cryst%symrec,cryst%symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk,ref_bz=ref_bz)
449 : else
450 784 : call identk(kibz,nkibz,nkbzX,nsym,timrev,cryst%symrec,cryst%symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk)
451 : end if
452 :
453 : ! TODO: Force the k-points to be in the first Brillouin zone.
454 : ! Now the GW tests seem to be OK, additional tests have to be done though.
455 806 : do_wrap=.FALSE.; if (PRESENT(wrap_1zone)) do_wrap=wrap_1zone
456 : !do_wrap=.TRUE.
457 :
458 357 : if (do_wrap) then ! Wrap the BZ points in the interval ]-1/2,1/2]
459 0 : do ik_bz=1,nkbz
460 0 : call wrap2_pmhalf(kbz(:,ik_bz),kbz_wrap,shift)
461 0 : kbz(:,ik_bz) = kbz_wrap
462 : end do
463 : end if
464 :
465 : ! ================================================================
466 : ! ==== Create data structure to store information on k-points ====
467 : ! ================================================================
468 : !
469 : ! Dimensions.
470 806 : Kmesh%nbz = nkbz ! Number of points in the full BZ
471 806 : Kmesh%nibz = nkibz ! Number of points in the IBZ
472 806 : Kmesh%nsym = nsym ! Number of operations
473 806 : Kmesh%timrev= timrev ! 2 if time-reversal is used, 1 otherwise
474 :
475 : ! Arrays.
476 10478 : Kmesh%gmet = Cryst%gmet
477 10478 : Kmesh%gprimd = Cryst%gprimd
478 :
479 2418 : ABI_MALLOC(Kmesh%bz ,(3,nkbz))
480 153116 : Kmesh%bz = kbz(:,1:nkbz ) ! Red. coordinates of points in full BZ.
481 2418 : ABI_MALLOC(Kmesh%ibz,(3,nkibz))
482 25092 : Kmesh%ibz = kibz(:,1:nkibz) ! Red. coordinates of points in IBZ.
483 :
484 2418 : ABI_MALLOC(Kmesh%tab ,(nkbz))
485 39488 : Kmesh%tab = ktab (1:nkbz) ! Index of the irred. point in the array IBZ.
486 1612 : ABI_MALLOC(Kmesh%tabi,(nkbz))
487 39488 : Kmesh%tabi = ktabi(1:nkbz) !-1 if time reversal must be used to obtain this point,
488 1612 : ABI_MALLOC(Kmesh%tabo,(nkbz))
489 39488 : Kmesh%tabo = ktabo(1:nkbz) ! Symm. operation that rotates k_IBZ onto \pm k_BZ
490 : ! (depending on tabi)
491 1612 : ABI_MALLOC(Kmesh%wt,(nkibz))
492 6676 : Kmesh%wt(:)= wtk(1:nkibz) ! Weight for each k_IBZ
493 :
494 4030 : ABI_MALLOC(Kmesh%rottbm1,(nkbz,timrev,nsym))
495 3224 : ABI_MALLOC(Kmesh%rottb ,(nkbz,timrev,nsym))
496 :
497 806 : do_hack = .FALSE.
498 806 : if (PRESENT(ref_bz) .and. PRESENT(break_symmetry)) then
499 22 : do_hack = break_symmetry
500 : end if
501 :
502 22 : if (do_hack) then
503 22 : ABI_WARNING("Hacking the rottb tables!")
504 1150 : do ik_bz=1,nkbz
505 163560 : Kmesh%rottbm1(ik_bz,:,:) = ik_bz
506 163582 : Kmesh%rottb (ik_bz,:,:) = ik_bz
507 : end do
508 : else
509 784 : call setup_k_rotation(nsym,timrev,cryst%symrec,nkbz,Kmesh%bz,Cryst%gmet,Kmesh%rottb,Kmesh%rottbm1)
510 : end if
511 :
512 : ! TODO umklp can be calculated inside setup_k_rotation.
513 2418 : ABI_MALLOC(Kmesh%umklp,(3,nkbz))
514 38682 : do ik_bz=1,nkbz
515 37876 : ik_ibz= Kmesh%tab (ik_bz)
516 37876 : isym = Kmesh%tabo(ik_bz)
517 37876 : itim = (3-Kmesh%tabi(ik_bz))/2
518 1250714 : Kmesh%umklp(:,ik_bz) = NINT( -Kmesh%bz(:,ik_bz) + (3-2*itim)*MATMUL(cryst%symrec(:,:,isym),Kmesh%ibz(:,ik_ibz)) )
519 : end do
520 :
521 2418 : ABI_MALLOC(Kmesh%tabp,(nkbz))
522 38682 : do ik_bz=1,nkbz
523 37876 : isym =Kmesh%tabo(ik_bz)
524 37876 : ik_ibz=Kmesh%tab (ik_bz)
525 946900 : rm1t=MATMUL(TRANSPOSE(cryst%symrec(:,:,isym)),cryst%tnons(:,isym))
526 152310 : Kmesh%tabp(ik_bz)=EXP(-(0.,1.)*two_pi*DOT_PRODUCT(kibz(:,ik_ibz),rm1t))
527 : end do
528 :
529 806 : ABI_FREE(kbz)
530 806 : ABI_FREE(wtk)
531 806 : ABI_FREE(ktab)
532 806 : ABI_FREE(ktabi)
533 806 : ABI_FREE(ktabo)
534 :
535 806 : end subroutine kmesh_init
536 : !!***
537 :
538 : !----------------------------------------------------------------------
539 :
540 : !!****f* m_bz_mesh/kmesh_free
541 : !! NAME
542 : !! kmesh_free
543 : !!
544 : !! FUNCTION
545 : !! Deallocate all dynamics entities present in a kmesh_t structure.
546 : !!
547 : !! SOURCE
548 :
549 806 : subroutine kmesh_free(Kmesh)
550 :
551 : !Arguments ------------------------------------
552 : class(kmesh_t),intent(inout) :: Kmesh
553 : ! *********************************************************************
554 :
555 : ! integer
556 806 : ABI_SFREE(Kmesh%rottb)
557 806 : ABI_SFREE(Kmesh%rottbm1)
558 806 : ABI_SFREE(Kmesh%tab)
559 806 : ABI_SFREE(Kmesh%tabi)
560 806 : ABI_SFREE(Kmesh%tabo)
561 806 : ABI_SFREE(Kmesh%umklp)
562 :
563 : ! real
564 806 : ABI_SFREE(Kmesh%ibz)
565 806 : ABI_SFREE(Kmesh%bz)
566 806 : ABI_SFREE(Kmesh%shift)
567 806 : ABI_SFREE(Kmesh%wt)
568 :
569 : ! complex
570 806 : ABI_SFREE(Kmesh%tabp)
571 :
572 806 : end subroutine kmesh_free
573 : !!***
574 :
575 : !----------------------------------------------------------------------
576 :
577 : !!****f* m_bz_mesh/kmesh_print
578 : !! NAME
579 : !! kmesh_print
580 : !!
581 : !! FUNCTION
582 : !! Print the content of a kmesh_t datatype
583 : !!
584 : !! INPUTS
585 : !! units: unit numbers
586 : !! [header]=optional header
587 : !! [prtvol]=verbosity level
588 : !!
589 : !! OUTPUT
590 : !! Only printing.
591 : !!
592 : !! SOURCE
593 :
594 1596 : subroutine kmesh_print(Kmesh, units, header, prtvol)
595 :
596 : !Arguments ------------------------------------
597 : !scalars
598 : class(kmesh_t),intent(in) :: Kmesh
599 : integer,intent(in) :: units(:)
600 : character(len=*),optional,intent(in) :: header
601 : integer,optional,intent(in) :: prtvol
602 :
603 : !Local variables-------------------------------
604 : !scalars
605 : integer,parameter :: nmaxk=50
606 : integer :: ii, ik, my_prtvol
607 : character(len=500) :: msg, fmt
608 : ! *************************************************************************
609 :
610 798 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
611 :
612 798 : msg=' ==== Info on the Kmesh% object ==== '
613 798 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
614 798 : call wrtout(units, msg)
615 :
616 : write(msg,'(a,i5,3a)')&
617 798 : ' Number of points in the irreducible wedge : ',Kmesh%nibz,ch10,&
618 1596 : ' Reduced coordinates and weights : ',ch10
619 798 : call wrtout(units, msg)
620 :
621 798 : write(fmt,*)'(1x,i5,a,2x,3es16.8,3x,f11.5)'
622 6609 : do ik=1,Kmesh%nibz ! Add tol8 for portability reasons.
623 23244 : write(msg,fmt) ik,') ',(Kmesh%ibz(ii,ik),ii=1,3),Kmesh%wt(ik)+tol8
624 6609 : call wrtout(units, msg)
625 : end do
626 :
627 828 : select case (Kmesh%timrev)
628 : case (1)
629 30 : write(msg,'(2a,i2,3a,i5,a)')ch10,&
630 30 : ' Together with ',Kmesh%nsym,' symmetry operations (time-reversal symmetry not used) ',ch10,&
631 60 : ' yields ',Kmesh%nbz,' points in the full Brillouin Zone.'
632 :
633 : case (2)
634 768 : write(msg,'(2a,i2,3a,i5,a)')ch10,&
635 768 : ' Together with ',Kmesh%nsym,' symmetry operations and time-reversal symmetry ',ch10,&
636 1536 : ' yields ',Kmesh%nbz,' points in the full Brillouin Zone.'
637 :
638 : case default
639 798 : ABI_BUG(sjoin('Wrong value for timrev:', itoa(Kmesh%timrev)))
640 : end select
641 :
642 798 : call wrtout(units, msg)
643 :
644 798 : if (my_prtvol > 0) then
645 104 : write(fmt,*)'(1x,i5,a,2x,3es16.8)'
646 4200 : do ik=1,Kmesh%nbz
647 4128 : if (my_prtvol==1 .and. ik>nmaxk) then
648 32 : call wrtout(units, ' prtvol=1, do not print more points.'); EXIT
649 : end if
650 16384 : write(msg,fmt)ik,') ',(Kmesh%bz(ii,ik),ii=1,3)
651 4168 : call wrtout(units, msg)
652 : end do
653 : end if
654 :
655 : ! Additional printing
656 104 : if (my_prtvol >= 10) then
657 0 : write(msg,'(2a)')ch10,&
658 0 : ' Full point -------> Irred point --> through: Symrec Time-Rev (1=No,-1=Yes) G0(1:3) '
659 0 : call wrtout(units, msg)
660 0 : write(fmt,*)'(2x,i5,2x,2(3(f7.4,2x)),i3,2x,i2,3(i3))'
661 0 : do ik=1,Kmesh%nbz
662 0 : write(msg,fmt) ik,Kmesh%bz(:,ik),Kmesh%ibz(:,Kmesh%tab(ik)),Kmesh%tabo(ik),Kmesh%tabi(ik),Kmesh%umklp(:,ik)
663 0 : call wrtout(units, msg)
664 : end do
665 : end if
666 :
667 798 : write(msg,'(a)')ch10
668 798 : call wrtout(units, msg)
669 :
670 798 : end subroutine kmesh_print
671 : !!***
672 :
673 : !----------------------------------------------------------------------
674 :
675 : !!****f* m_bz_mesh/setup_k_rotation
676 : !! NAME
677 : !! setup_k_rotation
678 : !!
679 : !! FUNCTION
680 : !! Set up tables giving the correspondence btw a k-point and its rotated image.
681 : !!
682 : !! INPUTS
683 : !! timrev=2 if time-reversal can be used, 1 otherwise.
684 : !! nsym=Number of symmetry operations
685 : !! symrec(3,3,nsym)=Symmetry operations in reciprocal space in reduced coordinates.
686 : !! nbz=Number of k-points
687 : !! kbz(3,nbz)=k-points in reduced coordinates.
688 : !! gmet(3,3)=Metric in reciprocal space.
689 : !!
690 : !! OUTPUT
691 : !! krottb(k,I,S)=Index of (IS) k in the array bz
692 : !! krottbm1(k,I,S)=Index of IS^{-1} k
693 : !!
694 : !! SOURCE
695 :
696 784 : subroutine setup_k_rotation(nsym, timrev, symrec, nbz, kbz, gmet, krottb, krottbm1)
697 :
698 : !Arguments ------------------------------------
699 : !scalars
700 : integer,intent(in) :: nbz,nsym,timrev
701 : !arrays
702 : integer,intent(in) :: symrec(3,3,nsym)
703 : integer,intent(out) :: krottb(nbz,timrev,nsym),krottbm1(nbz,timrev,nsym)
704 : real(dp),intent(in) :: kbz(3,nbz),gmet(3,3)
705 :
706 : !Local variables ------------------------------
707 : !scalars
708 : integer :: ik,ikp,isym,itim,nsh,ik_st !,sh_start
709 : real(dp),parameter :: KTOL=tol6
710 : real(dp) :: norm_old,norm !,norm_rot !norm_base
711 : logical :: found,isok
712 : character(len=500) :: msg
713 : !arrays
714 1568 : integer :: g0(3), iperm(nbz),shlim(nbz+1)
715 1568 : real(dp) :: shift(3),kbase(3),krot(3),knorm(nbz),kwrap(3),shlen(nbz+1)
716 : !************************************************************************
717 :
718 : DBG_ENTER("COLL")
719 :
720 : ! Sort the k-points according to their norm to speed up the search below.
721 37532 : do ik=1,nbz
722 146992 : call wrap2_pmhalf(kbz(:,ik), kwrap, shift)
723 36748 : knorm(ik) = normv(kwrap,gmet,"G")
724 37532 : iperm(ik)= ik
725 : end do
726 :
727 784 : call sort_dp(nbz,knorm,iperm,KTOL)
728 : !
729 : ! The index of the initial (sorted) k-point in each shell
730 784 : nsh=1; norm_old=knorm(1)
731 :
732 784 : shlim(1)=1; shlen(1) = norm_old
733 36748 : do ik_st=2,nbz
734 35964 : norm = knorm(ik_st)
735 36748 : if (ABS(norm-norm_old) > KTOL) then
736 5215 : norm_old = norm
737 5215 : nsh = nsh+1
738 5215 : shlim(nsh) = ik_st
739 5215 : shlen(nsh) = norm
740 : end if
741 : end do
742 784 : shlim(nsh+1)=nbz+1
743 784 : shlen(nsh+1)=HUGE(one)
744 : !
745 : ! === Set up k-rotation tables ===
746 : ! * Use spatial inversion instead of time reversal whenever possible.
747 : !call wrtout(std_out," Begin sorting ","COLL")
748 :
749 784 : isok=.TRUE.
750 37532 : do ik=1,nbz
751 146992 : kbase(:)=kbz(:,ik)
752 : !
753 110788 : do itim=1,timrev
754 3452356 : do isym=1,nsym
755 93585856 : krot(:)=(3-2*itim)*MATMUL(symrec(:,:,isym),kbase)
756 :
757 206529832 : found=.FALSE.
758 : #if 1
759 : ! Old code
760 206529832 : do ikp=1,nbz
761 206529832 : if (isamek(krot,kbz(:,ikp),g0)) then
762 3342352 : found=.TRUE.
763 3342352 : krottb (ik ,itim,isym)=ikp
764 3342352 : krottbm1(ikp,itim,isym)=ik
765 : EXIT
766 : end if
767 : end do
768 : #else
769 : ! Locate the shell index with bisection.
770 : call wrap2_pmhalf(krot(:),kwrap,shift)
771 : norm_rot = normv(kwrap,gmet,"G")
772 : sh_start = bisect(shlen(1:nsh+1),norm_rot)
773 :
774 : do ik_st=shlim(sh_start),nbz
775 : ikp = iperm(ik_st)
776 : if (isamek(krot,kbz(:,ikp),g0)) then
777 : found=.TRUE.
778 : krottb (ik ,itim,isym)=ikp
779 : krottbm1(ikp,itim,isym)=ik
780 : !write(std_out,*)ik_st,shlim(sh_start),nbz
781 : EXIT
782 : end if
783 : end do
784 : #endif
785 73256 : if (.not.found) then
786 0 : isok=.FALSE.
787 : !write(std_out,*)" norm_base,norm_rot ",norm_base,norm_rot
788 : !write(std_out,*)normv(kbase,gmet,"G"),normv(krot,gmet,"G")
789 : write(msg,'(2(a,i4),2x,2(3f12.6,2a),i3,a,i2)')&
790 0 : 'Initial k-point ',ik,'/',nbz,kbase(:),ch10,&
791 0 : 'Rotated k-point (not found) ',krot(:),ch10,&
792 0 : 'Through symmetry operation ',isym,' and itim ',itim
793 0 : ABI_ERROR(msg)
794 : end if
795 :
796 : end do
797 : end do
798 : end do
799 :
800 784 : if (.not.isok) then
801 0 : ABI_ERROR('k-mesh not closed')
802 : end if
803 :
804 : DBG_EXIT("COLL")
805 :
806 784 : end subroutine setup_k_rotation
807 : !!***
808 :
809 : !----------------------------------------------------------------------
810 :
811 : !!****f* m_bz_mesh/get_bz_item
812 : !! NAME
813 : !! get_bz_item
814 : !!
815 : !! FUNCTION
816 : !! Given the index of a point in the full BZ, this routine returns the index of the
817 : !! symmetric image in the IBZ, the index of the symmetry operation symrec needed,
818 : !! whether time-reversal has to be used.
819 : !! Optionally the non-symmorphic phase and the umklapp vector is returned.
820 : !!
821 : !! INPUTS
822 : !! ikbz=The index of the required point in the BZ
823 : !! Kmesh<kmesh_t>=Datatype gathering information on the k point sampling.
824 : !!
825 : !! OUTPUT
826 : !! kbz(3)=The k-point in the first BZ in reduced coordinated.
827 : !! isym=Index of the symrec symmetry required to rotate ik_ibz onto ik_bz.
828 : !! itim=2 is time-reversal has to be used, 1 otherwise
829 : !! ik_ibz=The index of the corresponding symmetric point in the IBZ.
830 : !! [ph_mkbzt]=The phase factor for non-symmorphic operations e^{-i 2 \pi k_IBZ \cdot R{^-1}t}=e{-i 2\pi k_BZ cdot t}
831 : !! [umklp(3)]=The umklapp G0 vector such as kbz + G0 = (IS) k_ibz, where kbz is in the BZ.
832 : !! [isirred]=.TRUE. if the k-point belongs to IBZ.
833 : !!
834 : !! SOURCE
835 :
836 751459 : subroutine get_bz_item(Kmesh, ik_bz, kbz, ik_ibz, isym, itim, ph_mkbzt, umklp, isirred)
837 :
838 : !Arguments ------------------------------------
839 : !scalars
840 : class(kmesh_t),intent(in) :: Kmesh
841 : integer,intent(in) :: ik_bz
842 : integer,intent(out) :: ik_ibz,isym,itim
843 : complex(dp),optional,intent(out) :: ph_mkbzt
844 : logical,optional,intent(out) :: isirred
845 : !arrays
846 : integer,optional,intent(out) :: umklp(3)
847 : real(dp),intent(out) :: kbz(3)
848 :
849 : !Local variables-------------------------------
850 : !scalars
851 : character(len=500) :: msg
852 : ! *********************************************************************
853 :
854 751459 : if (ik_bz>Kmesh%nbz.or.ik_bz<=0) then
855 0 : write(msg,'(a,2i3)')' Wrong value for ik_bz: ',ik_bz,Kmesh%nbz
856 0 : ABI_BUG(msg)
857 : end if
858 :
859 3005836 : kbz = Kmesh%bz(:,ik_bz)
860 751459 : ik_ibz = Kmesh%tab(ik_bz)
861 751459 : isym = Kmesh%tabo(ik_bz)
862 751459 : itim = (3-Kmesh%tabi(ik_bz))/2
863 :
864 751459 : if (PRESENT(ph_mkbzt)) ph_mkbzt=Kmesh%tabp(ik_bz)
865 1157479 : if (PRESENT(umklp)) umklp =Kmesh%umklp(:,ik_bz)
866 : ! Be careful here as we assume a particular ordering of symmetries.
867 2745506 : if (PRESENT(isirred)) isirred = (isym==1.and.itim==1.and.ALL(Kmesh%umklp(:,ik_bz)==(/0,0,0/)))
868 :
869 751459 : end subroutine get_bz_item
870 : !!***
871 :
872 : !----------------------------------------------------------------------
873 :
874 : !!****f* m_bz_mesh/get_IBZ_item
875 : !! NAME
876 : !! get_IBZ_item
877 : !!
878 : !! FUNCTION
879 : !! Report useful information on a k-point in the IBZ starting from its sequential index in %ibz.
880 : !!
881 : !! INPUTS
882 : !! ik_ibz=The index of the required point in the IBZ
883 : !! Kmesh<kmesh_t>=datatype gathering information on the k point sampling.
884 : !!
885 : !! OUTPUT
886 : !! kibz(3)=the k-point in reduced coordinated
887 : !! wtk=the weight
888 : !!
889 : !! TODO
890 : !! Add mapping ibz2bz, ibz2star
891 : !!
892 : !! SOURCE
893 :
894 308100 : subroutine get_IBZ_item(Kmesh,ik_ibz,kibz,wtk)
895 :
896 : !Arguments ------------------------------------
897 : !scalars
898 : class(kmesh_t),intent(in) :: Kmesh
899 : integer,intent(in) :: ik_ibz
900 : real(dp),intent(out) :: wtk
901 : !arrays
902 : real(dp),intent(out) :: kibz(3)
903 : ! *********************************************************************
904 :
905 308100 : if (ik_ibz>Kmesh%nibz.or.ik_ibz<=0) then
906 0 : ABI_BUG(sjoin('wrong value for ik_ibz: ',itoa(ik_ibz)))
907 : end if
908 :
909 1232400 : kibz=Kmesh%ibz(:,ik_ibz)
910 308100 : wtk =Kmesh%wt(ik_ibz)
911 :
912 308100 : end subroutine get_IBZ_item
913 : !!***
914 :
915 : !----------------------------------------------------------------------
916 :
917 : !!****f* m_bz_mesh/get_BZ_diff
918 : !! NAME
919 : !! get_BZ_diff
920 : !!
921 : !! FUNCTION
922 : !! Given two points k1 and k2 where k1 belongs to the BZ, check if the difference
923 : !! k1-k2 still belongs to the BZ reporting useful quantities
924 : !!
925 : !! INPUTS
926 : !! Kmesh<kmesh_t>=datatype gathering information on the k-mesh
927 : !! k1(3)=the first k-points (supposed to be in the BZ)
928 : !! k2(3)=the second point
929 : !!
930 : !! OUTPUT
931 : !! idiff_bz=the idex of k1-k2 in the BZ
932 : !! G0(3)=the umklapp G0 vector required to bring k1-k2 back to the BZ
933 : !! nfound= the number of points in the BZ that are equal to k1-k2 (should be 1 if everything is OK)
934 : !!
935 : !! SOURCE
936 :
937 67414 : subroutine get_BZ_diff(Kmesh,k1,k2,idiff_bz,g0,nfound)
938 :
939 : !Arguments ------------------------------------
940 : !scalars
941 : class(kmesh_t),intent(in) :: Kmesh
942 : integer,intent(out) :: idiff_bz,nfound
943 : !arrays
944 : integer,intent(out) :: g0(3)
945 : real(dp),intent(in) :: k1(3),k2(3)
946 :
947 : !Local variables-------------------------------
948 : !scalars
949 : integer :: ikp
950 : character(len=500) :: msg
951 : !arrays
952 : integer :: umklp(3)
953 : real(dp) :: kdiff(3),ktrial(3)
954 : ! *********************************************************************
955 :
956 67414 : if (.not.has_BZ_item(Kmesh,k1,ikp,umklp)) then
957 0 : write(msg,'(a,3f12.6)')' first point must be in BZ: ',k1
958 0 : ABI_ERROR(msg)
959 : end if
960 :
961 269656 : kdiff = k1-k2
962 67414 : nfound = 0
963 67414 : idiff_bz= 0
964 :
965 : ! === Find p such k1-k2=p+g0 where p in the BZ ===
966 11509906 : do ikp=1,Kmesh%nbz
967 45769968 : ktrial=Kmesh%bz(:,ikp)
968 11509906 : if (isamek(kdiff,ktrial,umklp)) then
969 67414 : idiff_bz=ikp
970 67414 : g0=umklp
971 67414 : nfound=nfound+1
972 : end if
973 : end do
974 :
975 : ! === Check if p has not found of found more than once ===
976 : ! * For extremely dense meshes, tol1q in defs_basis might be too large!
977 67414 : if (nfound/=1) then
978 0 : if (nfound==0) then
979 0 : ABI_WARNING(" k1-k2-G0 not found in BZ")
980 : else
981 0 : ABI_WARNING(sjoin(' Multiple k1-k2-G0 found in BZ, nfound= ', itoa(nfound)))
982 : end if
983 : write(msg,'(4a,3(a,3f12.6,a))') &
984 0 : ' k1 = ',k1 ,ch10,&
985 0 : ' k2 = ',k2 ,ch10,&
986 0 : ' k1-k2 = ',kdiff,ch10
987 0 : ABI_WARNING(msg)
988 : end if
989 :
990 67414 : end subroutine get_BZ_diff
991 : !!***
992 :
993 : !----------------------------------------------------------------------
994 :
995 : !!****f* m_bz_mesh/isamek
996 : !! NAME
997 : !! isamek
998 : !!
999 : !! FUNCTION
1000 : !! Test two k-points for equality.
1001 : !! Return .TRUE. is they are equal within a reciprocal lattice vector G0.
1002 : !!
1003 : !! INPUTS
1004 : !! k1(3),k2(3)=The two k points to be compared.
1005 : !!
1006 : !! OUTPUT
1007 : !! Return .TRUE. if they are the same within a RL vector,
1008 : !! .FALSE. if they are different.
1009 : !! G0(3)=if .TRUE. G0(3) is the reciprocal lattice vector such as k1=k2+G0
1010 : !!
1011 : !! SOURCE
1012 :
1013 465867233 : logical function isamek(k1, k2, g0)
1014 :
1015 : !Arguments ------------------------------------
1016 : !arrays
1017 : integer,intent(out) :: g0(3)
1018 : real(dp),intent(in) :: k1(3),k2(3)
1019 : ! *************************************************************************
1020 :
1021 1863468932 : isamek = isinteger(k1 - k2, TOL_KDIFF)
1022 :
1023 465867233 : if (isamek) then
1024 19561488 : g0 = NINT(k1-k2)
1025 : else
1026 1843907444 : g0 = HUGE(1)
1027 : end if
1028 :
1029 465867233 : end function isamek
1030 : !!***
1031 :
1032 : !----------------------------------------------------------------------
1033 :
1034 : !!****f* m_bz_mesh/isequalk
1035 : !! NAME
1036 : !! is_equalk
1037 : !!
1038 : !! FUNCTION
1039 : !! Return .TRUE. if two points are equal within a reciprocal lattice vector.
1040 : !!
1041 : !! INPUTS
1042 : !! q1(3),q2(3)=The two points to be compared for equivalence.
1043 : !!
1044 : !! OUTPUT
1045 : !!
1046 : !! SOURCE
1047 :
1048 1652964 : logical function isequalk(q1, q2)
1049 :
1050 : !Arguments ------------------------------------
1051 : real(dp),intent(in) :: q1(3),q2(3)
1052 :
1053 : !Local variables-------------------------------
1054 : integer :: g0(3)
1055 : ! *************************************************************************
1056 :
1057 60 : isequalk = isamek(q1,q2,g0)
1058 :
1059 60 : end function isequalk
1060 : !!***
1061 :
1062 : !----------------------------------------------------------------------
1063 :
1064 : !!****f* m_bz_mesh/has_BZ_item
1065 : !! NAME
1066 : !! has_BZ_item
1067 : !!
1068 : !! FUNCTION
1069 : !! check if item belongs to the BZ within a reciprocal lattice vector
1070 : !! and return the index number and the reciprocal vector g0.
1071 : !!
1072 : !! INPUTS
1073 : !! Kmesh<kmesh_t>=datatype gathering information on the k-mesh
1074 : !! item(3)=the k-point to be checked
1075 : !!
1076 : !! OUTPUT
1077 : !! .TRUE. if item is the BZ within a RL vector
1078 : !! ikbz=Index of the k-point in the Kmesh%bz array
1079 : !! g0(3)=Umklapp vector.
1080 : !!
1081 : !! FIXME
1082 : !! Switch to routine version. Due to side-effects the present implementation
1083 : !! might be source of bugs in logical statements
1084 : !!
1085 : !! SOURCE
1086 :
1087 74760 : logical function has_BZ_item(Kmesh, item, ikbz, g0)
1088 :
1089 : !Arguments ------------------------------------
1090 : !scalars
1091 : class(kmesh_t),intent(in) :: Kmesh
1092 : integer,intent(out) :: ikbz
1093 : !arrays
1094 : integer,intent(out) :: g0(3)
1095 : real(dp),intent(in) :: item(3)
1096 :
1097 : !Local variables-------------------------------
1098 : !scalars
1099 : integer :: ik_bz,yetfound
1100 : !arrays
1101 : integer :: g0_tmp(3)
1102 : ! *************************************************************************
1103 :
1104 74760 : has_BZ_item=.FALSE.; ikbz=0; g0=0; yetfound=0
1105 12857018 : do ik_bz=1,Kmesh%nbz
1106 12857018 : if (isamek(item, Kmesh%bz(:,ik_bz), g0_tmp)) then
1107 74760 : has_BZ_item=.TRUE.
1108 74760 : ikbz=ik_bz
1109 74760 : g0 = g0_tmp
1110 74760 : yetfound=yetfound+1
1111 : !EXIT
1112 : end if
1113 : end do
1114 :
1115 74760 : if (yetfound/=0 .and. yetfound/=1) then
1116 0 : ABI_ERROR('Multiple k-points found')
1117 : end if
1118 :
1119 74760 : end function has_BZ_item
1120 : !!***
1121 :
1122 : !----------------------------------------------------------------------
1123 :
1124 : !!****f* m_bz_mesh/has_IBZ_item
1125 : !! NAME
1126 : !! has_IBZ_item
1127 : !!
1128 : !! FUNCTION
1129 : !! Check if item belongs to the IBZ within a reciprocal lattice vector
1130 : !!
1131 : !! INPUTS
1132 : !! Kmesh<kmesh_t>=Datatype gathering information on the mesh in the BZ.
1133 : !! item(3)=the k-point to be checked
1134 : !!
1135 : !! OUTPUT
1136 : !! Return .TRUE. if item is the IBZ within a RL vector
1137 : !! ikibz=The index of the k-point in the IBZ.
1138 : !! g0(3)=The reciprocal lattice vector.
1139 : !!
1140 : !! SOURCE
1141 :
1142 517 : logical function has_IBZ_item(Kmesh, item, ikibz, g0)
1143 :
1144 : !Arguments ------------------------------------
1145 : !scalars
1146 : class(kmesh_t),intent(in) :: Kmesh
1147 : integer,intent(out) :: ikibz
1148 : !arrays
1149 : integer,intent(out) :: g0(3)
1150 : real(dp),intent(in) :: item(3)
1151 :
1152 : !Local variables-------------------------------
1153 : !scalars
1154 : integer :: ik_ibz,yetfound
1155 : !character(len=500) :: msg
1156 : !arrays
1157 : integer :: g0_tmp(3)
1158 : ! *************************************************************************
1159 :
1160 517 : has_IBZ_item=.FALSE.; ikibz=0; g0=0; yetfound=0
1161 3994 : do ik_ibz=1,Kmesh%nibz
1162 3994 : if (isamek(item, Kmesh%ibz(:,ik_ibz), g0_tmp)) then
1163 517 : has_IBZ_item = .TRUE.
1164 517 : ikibz=ik_ibz
1165 517 : g0 = g0_tmp
1166 517 : yetfound=yetfound+1
1167 : !EXIT
1168 : end if
1169 : end do
1170 :
1171 517 : if (yetfound /=0 .and. yetfound /= 1) then
1172 0 : ABI_BUG("multiple k-points found")
1173 : end if
1174 :
1175 517 : end function has_IBZ_item
1176 : !!***
1177 :
1178 : !----------------------------------------------------------------------
1179 :
1180 : !!****f* m_bz_mesh/bz_mesh_isirred
1181 : !! NAME
1182 : !! bz_mesh_isirred
1183 : !!
1184 : !! FUNCTION
1185 : !! bz_mesh_isirred
1186 : !!
1187 : !! INPUTS
1188 : !! ik_bz=Index of the k-point in the BZ.
1189 : !!
1190 : !! OUTPUT
1191 : !! Return TRUE. if the k-point is in the IBZ (NB: a non-zero umklapp is not allowed)
1192 : !!
1193 : !! SOURCE
1194 :
1195 0 : pure logical function bz_mesh_isirred(Kmesh, ik_bz)
1196 :
1197 : !Arguments ------------------------------------
1198 : !scalars
1199 : class(kmesh_t),intent(in) :: Kmesh
1200 : integer,intent(in) :: ik_bz
1201 :
1202 : !Local variables-------------------------------
1203 : !scalars
1204 : integer :: isym,itim
1205 : ! *********************************************************************
1206 :
1207 0 : isym = Kmesh%tabo(ik_bz)
1208 0 : itim = (3-Kmesh%tabi(ik_bz))/2
1209 :
1210 : ! Be careful here as we assume a particular ordering of symmetries.
1211 0 : bz_mesh_isirred = (isym==1 .and. itim==1 .and. ALL(Kmesh%umklp(:,ik_bz) == [0,0,0]))
1212 :
1213 0 : end function bz_mesh_isirred
1214 : !!***
1215 :
1216 : !----------------------------------------------------------------------
1217 :
1218 : !!****f* m_bz_mesh/make_mesh
1219 : !! NAME
1220 : !! make_mesh
1221 : !!
1222 : !! FUNCTION
1223 : !! Initialize the kmesh_t starting from kptrlatt and shiftk
1224 : !!
1225 : !! INPUTS
1226 : !! Cryst<crystal_t>=Info on the crystalline structure.
1227 : !! nshiftk=Number of shifts for the mesh.
1228 : !! kptrlatt(3,3)= Coordinates of three vectors in real space, expressed in reduced coordinates.
1229 : !! They define a super-lattice in real space. The k point lattice is the reciprocal of
1230 : !! this super-lattice, eventually shifted by shift.
1231 : !! shiftk(3,nshiftk)=Shifts for the k-mesh.
1232 : !! [vacuum(3)]=For each direction, 0 if no vacuum, 1 if vacuum
1233 : !!
1234 : !! OUTPUT
1235 : !! Kmesh<kmesh_t>=Object gathering info on the sampling of the Brillouin zone.
1236 : !!
1237 : !! SOURCE
1238 :
1239 22 : subroutine make_mesh(Kmesh, Cryst, kptopt, kptrlatt, nshiftk, shiftk,&
1240 : vacuum, break_symmetry) ! Optional
1241 :
1242 : !Arguments -------------------------------
1243 : !scalars
1244 : type(kmesh_t),intent(inout) :: Kmesh
1245 : integer,intent(in) :: nshiftk,kptopt
1246 : logical,optional,intent(in) :: break_symmetry
1247 : type(crystal_t),intent(in) :: Cryst
1248 : !arrays
1249 : integer,intent(inout) :: kptrlatt(3,3)
1250 : integer,optional,intent(in) :: vacuum(3)
1251 : real(dp),intent(in) :: shiftk(3,nshiftk)
1252 :
1253 : !Local variables -------------------------
1254 : !scalars
1255 : integer,parameter :: chksymbreak0=0
1256 : integer :: iscf,nkbz,nkibz,nkpt_computed,my_nshiftk
1257 : real(dp) :: kptrlen
1258 : logical :: my_break_symmetry
1259 : !arrays
1260 : integer :: my_vacuum(3)
1261 22 : real(dp),allocatable :: kibz(:,:),wtk(:),my_shiftk(:,:),ref_kbz(:,:)
1262 : ! *************************************************************************
1263 :
1264 : DBG_ENTER("COLL")
1265 :
1266 22 : if (ALL(kptopt /= [1,3])) then
1267 0 : ABI_WARNING(sjoin("Not allowed value for kptopt: ", itoa(kptopt)))
1268 : end if
1269 : !
1270 : ! ======================================================================
1271 : ! ==== First call to getkgrid to obtain nkibz as well as the BZ set ====
1272 : ! ======================================================================
1273 22 : iscf=7 ! use for the Weights in NSCF calculation. check it more carefully.
1274 22 : nkibz=0 ! Compute number of k-points in the BZ and IBZ
1275 :
1276 22 : my_vacuum = [0,0,0]; if (PRESENT(vacuum)) my_vacuum=vacuum
1277 :
1278 22 : my_nshiftk = nshiftk
1279 22 : ABI_CHECK(my_nshiftk>0.and.my_nshiftk<=MAX_NSHIFTK, sjoin("Wrong nshiftk must be between 1 and ", itoa(MAX_NSHIFTK)))
1280 22 : ABI_MALLOC(my_shiftk, (3, MAX_NSHIFTK))
1281 18590 : my_shiftk=zero; my_shiftk(:,1:nshiftk) = shiftk(:,:)
1282 :
1283 : !write(std_out,*)" In make_mesh"; write(std_out,*)" kptopt ",kptopt," kptrlatt ",kptrlatt; !write(std_out,*)" nshiftk ",nshiftk," shiftk ",shiftk
1284 :
1285 22 : ABI_MALLOC(kibz,(3,nkibz))
1286 22 : ABI_MALLOC(wtk,(nkibz))
1287 :
1288 : call getkgrid(chksymbreak0,0,iscf,kibz,kptopt,kptrlatt,kptrlen,Cryst%nsym,0,nkibz,my_nshiftk,&
1289 22 : Cryst%nsym,Cryst%rprimd,my_shiftk,Cryst%symafm,Cryst%symrel,my_vacuum,wtk,fullbz=ref_kbz)
1290 :
1291 22 : nkbz = SIZE(ref_kbz,DIM=2)
1292 :
1293 22 : ABI_FREE(kibz)
1294 22 : ABI_FREE(wtk)
1295 :
1296 : !write(std_out,*)" after getkgrid1: nkbz = ",nkbz," nkibz=",nkibz; write(std_out,*)" ref_kbz = ",ref_kbz
1297 :
1298 : ! ==============================================================
1299 : ! ==== Recall getkgrid to get kibz(3,nkibz) and wtk(nkibz) =====
1300 : ! ==============================================================
1301 :
1302 66 : ABI_MALLOC(kibz,(3,nkibz))
1303 66 : ABI_MALLOC(wtk,(nkibz))
1304 :
1305 : call getkgrid(chksymbreak0,0,iscf,kibz,kptopt,kptrlatt,kptrlen,Cryst%nsym,nkibz,nkpt_computed,my_nshiftk,&
1306 22 : Cryst%nsym,Cryst%rprimd,my_shiftk,Cryst%symafm,Cryst%symrel,my_vacuum,wtk)
1307 :
1308 : ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ.
1309 22 : Kmesh%nshift = my_nshiftk
1310 286 : Kmesh%kptrlatt = kptrlatt
1311 :
1312 : ! Call the main creation method to get the tables tabo, tabi, tabp, umklp...
1313 : ! kmesh_init will reconstruct the BZ from kibz but pruning the k-points not in ref_bz
1314 : ! TODO: solve problem with timrev
1315 22 : my_break_symmetry=.FALSE.; if (PRESENT(break_symmetry)) my_break_symmetry=break_symmetry
1316 22 : call Kmesh%init(Cryst, nkibz, kibz, kptopt, ref_bz=ref_kbz, break_symmetry=my_break_symmetry)
1317 :
1318 22 : ABI_FREE(ref_kbz)
1319 22 : ABI_FREE(kibz)
1320 22 : ABI_FREE(wtk)
1321 :
1322 66 : ABI_MALLOC(Kmesh%shift,(3,my_nshiftk))
1323 132 : Kmesh%shift=my_shiftk(:,1:my_nshiftk)
1324 : ! Init Kmesh is breaking nshiftk
1325 22 : Kmesh%nshift=my_nshiftk
1326 :
1327 22 : ABI_FREE(my_shiftk)
1328 :
1329 : DBG_EXIT("COLL")
1330 :
1331 22 : end subroutine make_mesh
1332 : !!***
1333 :
1334 : !----------------------------------------------------------------------
1335 :
1336 : !!****f* m_bz_mesh/identk
1337 : !! NAME
1338 : !! identk
1339 : !!
1340 : !! FUNCTION
1341 : !! Identify k-points in the whole BZ starting from the IBZ.
1342 : !! Generate also symmetry tables relating the BZ to the IBZ.
1343 : !!
1344 : !! INPUTS
1345 : !! kibz(3,nkibz)=Coordinates of k-points in the IBZ.
1346 : !! nkibz=Number of k points in IBZ.
1347 : !! nkbzmx=Maximum number of k points in BZ.
1348 : !! nsym=Number of symmetry operations.
1349 : !! timrev=2 if time reversal symmetry can be used; 1 otherwise.
1350 : !! symrec(3,3,nsym)=Symmetry operation matrices in reciprocal space.
1351 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations.
1352 : !! [ref_bz(:,:)]= Reference set of points in the full Brillouin zone.
1353 : !!
1354 : !! OUTPUT
1355 : !! kbz(3,nkbzmx)= k-points in whole BZ
1356 : !! ktab(nkbzmx)= table giving for each k-point in the BZ (array kbz),
1357 : !! the corresponding irreducible point in the array (kibz)
1358 : !! k_BZ= (IS) kIBZ where S is one of the symrec operations and I is the inversion or the identity
1359 : !! where k_BZ = (IS) k_IBZ and S = \transpose R^{-1}
1360 : !! ktabi(nkbzmx)= for each k-point in the BZ defines whether inversion has to be
1361 : !! considered in the relation k_BZ=(IS) k_IBZ (1 => only S; -1 => -S)
1362 : !! ktabo(nkbzmx)= the symmetry operation S that takes k_IBZ to each k_BZ
1363 : !! nkbz= no. of k-points in the whole BZ
1364 : !! wtk(nkibz)= weight for each k-point in IBZ for symmetric quantities:
1365 : !! no. of distinct ks in whole BZ/(timrev*nsym)
1366 : !!
1367 : !! NOTES
1368 : !! The logic of the routine relies on the assumption that kibz really represent an irreducible set.
1369 : !! If symmetrical points are present in the input list, indeed, some the output weights will turn out to be zero.
1370 : !! An initial check is done at the beginning of the routine to trap this possible error.
1371 : !!
1372 : !! SOURCE
1373 :
1374 806 : subroutine identk(kibz, nkibz, nkbzmx, nsym, timrev, symrec, symafm, kbz, ktab, ktabi, ktabo, nkbz, wtk, ref_bz)
1375 :
1376 : !Arguments ------------------------------------
1377 : !scalars
1378 : integer,intent(in) :: nkbzmx,nkibz,nsym,timrev
1379 : integer,intent(out) :: nkbz
1380 : !arrays
1381 : integer,intent(in) :: symafm(nsym),symrec(3,3,nsym)
1382 : integer,intent(out) :: ktab(nkbzmx),ktabi(nkbzmx),ktabo(nkbzmx)
1383 : real(dp),intent(in) :: kibz(3,nkibz)
1384 : real(dp),intent(out) :: kbz(3,nkbzmx),wtk(nkibz)
1385 : real(dp),optional,intent(in) :: ref_bz(:,:)
1386 :
1387 : !Local variables ------------------------------
1388 : !scalars
1389 : integer :: ik1,ik2,ikbz,ikibz,iold,isym,itim
1390 : integer :: ikref,nkref,isym_swp,itim_swp,ikibz_swp
1391 : logical :: is_irred_set, found,ltest
1392 : character(len=500) :: msg
1393 : !arrays
1394 : integer :: g0(3)
1395 : real(dp) :: knew(3),k1(3),k2(3),kref(3),kbz_swp(3)
1396 : ! *************************************************************************
1397 :
1398 : DBG_ENTER("COLL")
1399 : !
1400 : ! === Check whether kibz really forms an irreducible set ===
1401 806 : is_irred_set=.TRUE.
1402 5870 : do ik1=1,nkibz-1
1403 20256 : k1=kibz(:,ik1)
1404 58710 : do ik2=ik1+1,nkibz
1405 211360 : k2=kibz(:,ik2)
1406 :
1407 163492 : do itim=1,timrev
1408 5140780 : do isym=1,nsym
1409 4982352 : if (symafm(isym)==-1) CYCLE
1410 138683328 : knew = (3-2*itim) * MATMUL(symrec(:,:,isym),k2)
1411 5058564 : if (isamek(k1,knew,g0)) then
1412 0 : is_irred_set=.FALSE.
1413 0 : write(msg,'(2(a,3f8.4),2(a,i0))')' k1 = ',k1,' is symmetrical of k2 = ',k2,' through sym = ',isym,' itim = ',itim
1414 0 : ABI_WARNING(msg)
1415 : end if
1416 : end do
1417 : end do
1418 :
1419 : end do
1420 : end do
1421 :
1422 : !call klist_isirred(nkibz,kibz,Cryst,nimg)
1423 :
1424 806 : if (.not.is_irred_set) then
1425 0 : ABI_WARNING("Input array kibz does not constitute an irreducible set.")
1426 : end if
1427 :
1428 : ! === Loop over k-points in IBZ ===
1429 : ! Start with zero no. of k-points found.
1430 806 : nkbz = 0
1431 6676 : do ikibz=1,nkibz
1432 5870 : wtk(ikibz) = zero
1433 :
1434 : ! === Loop over time-reversal I and symmetry operations S ===
1435 : ! Use spatial inversion instead of time reversal whenever possible.
1436 18348 : do itim=1,timrev
1437 539342 : do isym=1,nsym
1438 521800 : if (symafm(isym)==-1) CYCLE
1439 : !
1440 : ! * Form IS k
1441 14384608 : knew=(3-2*itim)*MATMUL(symrec(:,:,isym),kibz(:,ikibz))
1442 : !
1443 : ! * Check whether it has already been found (to within a RL vector).
1444 513736 : iold=0
1445 179758108 : do ikbz=1,nkbz
1446 179758108 : if (isamek(knew,kbz(:,ikbz),g0)) then
1447 513736 : iold=iold+1
1448 : exit
1449 : end if
1450 : end do
1451 : !
1452 : ! If not yet found add to kbz and increase the weight.
1453 11672 : if (iold==0) then
1454 90892 : nkbz=nkbz+1
1455 90892 : wtk(ikibz)=wtk(ikibz)+one
1456 90892 : if (nkbz>nkbzmx) then
1457 0 : ABI_BUG(sjoin('nkbzmx too small, nkbzmx = ',itoa(nkbzmx),', increase nkbzmx !'))
1458 : end if
1459 363568 : kbz(:,nkbz) = knew(:)
1460 90892 : ktab (nkbz) = ikibz
1461 90892 : ktabo(nkbz) = isym
1462 90892 : ktabi(nkbz) = 3-2*itim
1463 : end if
1464 : !
1465 : end do
1466 : end do
1467 :
1468 : end do !ikibz
1469 :
1470 806 : if (PRESENT(ref_bz)) then
1471 22 : call wrtout(std_out," Pruning the k-points not in ref_bz then reordering tables","COLL")
1472 :
1473 22 : nkref = SIZE(ref_bz,DIM=2)
1474 22 : ltest = (nkref>=nkbz.and.nkref<=nkbzmx)
1475 : if (.not.ltest) then
1476 22 : write(msg,'(3(a,i0))')" Wrong value for nkref: nkref= ",nkref," nkbz= ",nkbz," nkbzmx =",nkbzmx
1477 22 : ABI_WARNING(msg)
1478 : end if
1479 :
1480 1150 : do ikref=1,nkref
1481 4512 : kref = ref_bz(:,ikref)
1482 1128 : found=.FALSE.
1483 :
1484 1652904 : do ikbz=1,nkbz ! Loop on the set of BZ points found above.
1485 1652904 : if (isequalk(kref,kbz(:,ikbz))) then ! Swap indices.
1486 4512 : kbz_swp = kbz(:,ikref)
1487 1128 : ikibz_swp = ktab (ikref)
1488 1128 : isym_swp = ktabo(ikref)
1489 1128 : itim_swp = ktabi(ikref)
1490 :
1491 4512 : kbz(:,ikref) = kref
1492 1128 : ktab (ikref) = ktab (ikbz)
1493 1128 : ktabo(ikref) = ktabo(ikbz)
1494 1128 : ktabi(ikref) = ktabi(ikbz)
1495 :
1496 4512 : kbz(:,ikbz) = kbz_swp
1497 1128 : ktab (ikbz) = ikibz_swp
1498 1128 : ktabo(ikbz) = isym_swp
1499 1128 : ktabi(ikbz) = itim_swp
1500 :
1501 1128 : found=.TRUE.; EXIT
1502 : end if
1503 : end do
1504 :
1505 22 : if (.not.found) then
1506 0 : write(msg,'(a,3es16.8)')" One of the k-point in ref_bz is not a symmetrical image of the IBZ: ",kref
1507 0 : ABI_ERROR(msg)
1508 : end if
1509 : end do
1510 : !
1511 : ! Change nkbz to nkref, then get the new weights.
1512 1150 : nkbz=nkref; wtk=zero
1513 1150 : do ikref=1,nkref
1514 1128 : ikibz = ktab(ikref)
1515 1150 : wtk(ikibz) = wtk(ikibz) + 1
1516 : end do
1517 : end if ! PRESENT(ref_bz)
1518 : !
1519 : ! * Weights are normalized to 1.
1520 12546 : wtk = wtk/SUM(wtk)
1521 :
1522 : DBG_EXIT("COLL")
1523 :
1524 806 : end subroutine identk
1525 : !!***
1526 :
1527 : !----------------------------------------------------------------------
1528 :
1529 : !!****f* m_bz_mesh/get_ng0sh
1530 : !! NAME
1531 : !! get_ng0sh
1532 : !!
1533 : !! FUNCTION
1534 : !! Given two lists of k-points, kbz1 and kbz2, calculate any possible difference k1-k2.
1535 : !! For each difference, find the umklapp g0 vector and the point k3 in the array kfold
1536 : !! such as k1-k2 = k3 + G0.
1537 : !! The optimal value of G0 shells is returned, namely the smallest box around Gamma
1538 : !! which suffices to treat all possible umklapp processes.
1539 : !! The search algorithm uses bisection to process to scale in nk1*nk2*log(nkfold)
1540 : !!
1541 : !! INPUTS
1542 : !! nk1, nk2=Number of points in the arrays kbz1, kbz2.
1543 : !! kbz1(3,nk1)=Reduced coordinates of the first set of points.
1544 : !! kbz2(3,nk2)=Reduced coordinates of the second set of points.
1545 : !! nkfold=Number of points in the array kfold.
1546 : !! kfold(3,nkfol)=Reduced coordinated of the points in the BZ.
1547 : !! tolq0=Tolerance below which a q-point is treated as zero.
1548 : !!
1549 : !! OUTPUT
1550 : !! opt_ng0(3)=Minimal reduced components of the G0 vectors to account for umklapps.
1551 : !!
1552 : !! SOURCE
1553 :
1554 378 : subroutine get_ng0sh(nk1,kbz1,nk2,kbz2,nkfold,kfold,tolq0,opt_ng0)
1555 :
1556 : !Arguments ------------------------------------
1557 : !scalars
1558 : integer,intent(in) :: nk1,nk2,nkfold
1559 : real(dp),intent(in) :: tolq0
1560 : !arrays
1561 : integer,intent(out) :: opt_ng0(3)
1562 : real(dp),intent(in) :: kbz1(3,nk1),kbz2(3,nk2),kfold(3,nkfold)
1563 :
1564 : !Local variables-------------------------------
1565 : !scalars
1566 : integer :: i1,i2,ikf,ind,factor
1567 : real(dp) :: normdiff,tempnorm,smallestlen
1568 : logical :: found
1569 : character(len=500) :: msg
1570 : !arrays
1571 756 : integer :: roundk(3),kbigdiff(3),kbigfold(3,nkfold),iperm(nkfold)
1572 756 : real(dp) :: k1mk2(3),ksmalldiff(3),norm(nkfold),ksmallfold(3,nkfold)
1573 : !************************************************************************
1574 :
1575 : ! Compute smallest length of one component
1576 : ! To get a sufficiently large factor to order vectors
1577 378 : smallestlen = one
1578 :
1579 : ! Compute integer part and fractional part, [0,1[ of kfold
1580 18580 : do ikf = 1,nkfold
1581 72808 : kbigfold(:,ikf) = FLOOR(kfold(:,ikf)+tol7)
1582 72808 : ksmallfold(:,ikf) = kfold(:,ikf)-kbigfold(:,ikf)
1583 :
1584 18202 : if (ABS(ksmallfold(1,ikf)) > tol7) smallestlen = MIN(smallestlen, ABS(ksmallfold(1,ikf)))
1585 18202 : if (ABS(ksmallfold(2,ikf)) > tol7) smallestlen = MIN(smallestlen, ABS(ksmallfold(2,ikf)))
1586 18580 : if (ABS(ksmallfold(3,ikf)) > tol7) smallestlen = MIN(smallestlen, ABS(ksmallfold(3,ikf)))
1587 :
1588 : end do
1589 :
1590 : ! WARNING ! This could not be sufficient if tested k1 - k2 has lower
1591 : ! components than smallestlen. The factor 10 is giving us a security margin.
1592 378 : factor = 10*(int(one/smallestlen)+1)
1593 :
1594 : ! Loop again over kfold vectors to give each term its norm
1595 18580 : do ikf=1, nkfold
1596 18202 : iperm(ikf) = ikf
1597 :
1598 : ! Computing a sort of norm with order of components (used for ordering)
1599 18202 : call getkptnorm_bycomponent(ksmallfold(:,ikf),factor,tempnorm)
1600 :
1601 18580 : norm(ikf) = tempnorm
1602 : end do
1603 :
1604 : ! Sorting list of kfold vectors
1605 378 : call sort_dp(nkfold,norm,iperm,tol14)
1606 :
1607 : ! Loop over all k1 - k2
1608 378 : opt_ng0(:)=0
1609 12578 : do i2=1,nk2
1610 : ! This is used in case of screening calculation.
1611 : ! If q is small treat it as zero. In this case, indeed,
1612 : ! we use q=0 to calculate the oscillator matrix elements.
1613 12200 : if (is_zero(kbz2(:,i2),tolq0)) CYCLE
1614 161346 : do i1=1,nk1
1615 148916 : found=.FALSE.
1616 :
1617 : ! Separating in integer part and fractionary part
1618 595664 : k1mk2(:) = kbz1(:,i1)-kbz2(:,i2)
1619 : ! Adding small tol to prevent 1 being in fractionary part
1620 595664 : kbigdiff(:) = FLOOR(k1mk2(:)+tol7)
1621 595664 : ksmalldiff(:) = k1mk2(:)-kbigdiff(:)
1622 :
1623 148916 : call getkptnorm_bycomponent(ksmalldiff(:),factor,normdiff)
1624 :
1625 : ! Try to find the right smallkfold, corresponding to ksmalldiff
1626 148916 : ind = bisect(norm,normdiff)
1627 :
1628 148916 : if (ind > 0) then
1629 148590 : if(ABS(norm(ind) - normdiff) < TOL_KDIFF) then
1630 148916 : found = .TRUE.
1631 : end if
1632 : end if
1633 148916 : if(ind < nkfold) then
1634 148388 : if(ABS(norm(ind+1) - normdiff) < TOL_KDIFF) then
1635 : found = .TRUE.
1636 : ind = ind + 1
1637 : end if
1638 : end if
1639 :
1640 251522 : if (.not. found) then
1641 : write(msg,'(a,2(2a,i4,3es16.8),a)')&
1642 0 : 'Not able to found umklapp G0 vector such as k1-k2 = kf+G0',ch10,&
1643 0 : 'point1 = ',i1,kbz1(:,i1),ch10,&
1644 0 : 'point2 = ',i2,kbz2(:,i2),ch10
1645 0 : ABI_ERROR(msg)
1646 : else
1647 : ! We have found one k, extracting the max g0
1648 595664 : roundk(:) = ABS(kbigdiff(:) - kbigfold(:,iperm(ind)))
1649 148916 : opt_ng0(1) = MAX(opt_ng0(1),roundk(1))
1650 148916 : opt_ng0(2) = MAX(opt_ng0(2),roundk(2))
1651 148916 : opt_ng0(3) = MAX(opt_ng0(3),roundk(3))
1652 : end if
1653 : end do
1654 : end do
1655 :
1656 378 : end subroutine get_ng0sh
1657 : !!***
1658 :
1659 : !----------------------------------------------------------------------
1660 :
1661 : !!****f* m_bz_mesh/getkptnorm_bycomponent
1662 : !! NAME
1663 : !! getkptnorm_bycomponent
1664 : !!
1665 : !! FUNCTION
1666 : !! get the norm of one vector, by order of its components
1667 : !!
1668 : !! INPUTS
1669 : !! vect(3) = vector which should be evaluated
1670 : !! factor = term which multiplies the components
1671 : !! WARNING ! should be large enough to get unique order
1672 : !!
1673 : !! OUTPUT
1674 : !! norm = value of the norm
1675 : !!
1676 : !! SOURCE
1677 :
1678 167118 : subroutine getkptnorm_bycomponent(vect,factor,norm)
1679 :
1680 : !Arguments ------------------------------------
1681 : !scalars
1682 : integer,intent(in) :: factor
1683 : real(dp),intent(out):: norm
1684 : !arrays
1685 : real(dp),intent(in) :: vect(3)
1686 :
1687 : !Local variables-------------------------------
1688 : character(len=500) :: msg
1689 : ! *************************************************************************
1690 :
1691 : ! Checking the factor is large enough (skipping zero components, since in this case the product will be 0)
1692 668472 : if (ANY(vect(:)*factor < 1.0 .and. vect(:) > tol7)) then
1693 0 : write(msg,'(a,a,a,a,a,a,a,a)') ' Not able to give unique norm to order vectors',ch10,&
1694 0 : 'This is likely related to a truncation error for a k-point in the input file',ch10,&
1695 0 : 'Always prefer fractional numbers in the input file instead of truncated ones',ch10,&
1696 0 : '(e.g. 1/6 instead of 0.166666667)',ch10
1697 0 : ABI_ERROR(msg)
1698 : end if
1699 :
1700 167118 : norm = (vect(1)*factor+vect(2))*factor+vect(3)
1701 :
1702 167118 : end subroutine getkptnorm_bycomponent
1703 : !!***
1704 :
1705 : !----------------------------------------------------------------------
1706 :
1707 : !!****f* m_bz_mesh/make_path
1708 : !! NAME
1709 : !! make_path
1710 : !!
1711 : !! FUNCTION
1712 : !! Generate a normalized path given the extrema.
1713 : !! See also kpath_t and kpath_init (recommended API).
1714 : !!
1715 : !! INPUTS
1716 : !! nbounds=Number of extrema defining the path.
1717 : !! bounds(3,nbounds)=The points defining the path in reduced coordinates.
1718 : !! met(3,3)=Metric matrix.
1719 : !! space='R' for real space, G for reciprocal space.
1720 : !! ndivsm=Number of divisions to be used for the smallest segment.
1721 : !! [unit]=Fortran unit for formatted output. Default: dev_null
1722 : !!
1723 : !! OUTPUT
1724 : !! npts=Total number of points in the normalized circuit.
1725 : !! ndivs(nbounds-1)=Number of division for each segment
1726 : !! path: allocated inside the routine. When the subroutine returns, path(3,npts) will
1727 : !! contain the path in reduced coordinates.
1728 : !!
1729 : !! SOURCE
1730 :
1731 80 : subroutine make_path(nbounds, bounds, met, space, ndivsm, ndivs, npts, path, unit)
1732 :
1733 : !Arguments ------------------------------------
1734 : !scalars
1735 : integer,intent(in) :: nbounds,ndivsm
1736 : integer,optional,intent(in) :: unit
1737 : integer,intent(out) :: npts
1738 : character(len=1),intent(in) :: space
1739 : !arrays
1740 : integer,intent(out) :: ndivs(nbounds-1)
1741 : real(dp),intent(in) :: bounds(3,nbounds),met(3,3)
1742 : real(dp),allocatable,intent(out) :: path(:,:)
1743 :
1744 : !Local variables-------------------------------
1745 : !scalars
1746 : integer,parameter :: prtvol=0
1747 : integer :: idx,ii,jp,ount
1748 : real(dp) :: nfact
1749 : character(len=500) :: msg
1750 : !arrays
1751 160 : real(dp) :: diff(3),lng(nbounds-1)
1752 : ! *************************************************************************
1753 :
1754 80 : ABI_CHECK(ndivsm > 0, sjoin('ndivsm', itoa(ndivsm)))
1755 :
1756 80 : ount = dev_null; if (present(unit)) ount = unit
1757 :
1758 418 : do ii=1,nbounds-1
1759 1352 : diff(:)=bounds(:,ii+1)-bounds(:,ii)
1760 418 : lng(ii) = normv(diff,met,space)
1761 : end do
1762 :
1763 : ! Avoid division by zero if any k(:,i+1)=k(:,i).
1764 498 : nfact=MINVAL(lng)
1765 80 : if (ABS(nfact)<tol6) then
1766 : write(msg,'(3a)')&
1767 0 : 'Found two equivalent consecutive points in the path ',ch10,&
1768 0 : 'This is not allowed, modify the path in your input file'
1769 0 : ABI_ERROR(msg)
1770 : end if
1771 :
1772 80 : nfact=nfact/ndivsm
1773 418 : ndivs(:)=NINT(lng(:)/nfact)
1774 418 : npts=SUM(ndivs)+1 !1 for the first point
1775 :
1776 80 : write(msg,'(2a,i0,2a)')ch10,&
1777 80 : ' Total number of points in the path: ',npts,ch10,&
1778 160 : ' Number of divisions for each segment of the normalized path: '
1779 80 : call wrtout(ount,msg)
1780 :
1781 418 : do ii=1,nbounds-1
1782 338 : write(msg,'(2(3f8.5,a),i0,a)')bounds(:,ii),' ==> ',bounds(:,ii+1),' ( ndivs : ',ndivs(ii),' )'
1783 418 : call wrtout(ount,msg)
1784 : end do
1785 80 : call wrtout(ount,ch10)
1786 :
1787 : ! Allocate and construct the path.
1788 240 : ABI_MALLOC(path,(3,npts))
1789 :
1790 : if (prtvol > 0) call wrtout(ount,' Normalized Path: ')
1791 80 : idx=0
1792 418 : do ii=1,nbounds-1
1793 9554 : do jp=1,ndivs(ii)
1794 9136 : idx=idx+1
1795 36544 : path(:,idx)=bounds(:,ii)+(jp-1)*(bounds(:,ii+1)-bounds(:,ii))/ndivs(ii)
1796 338 : if (prtvol > 0) then
1797 : write(msg,'(i4,4x,3(f8.5,1x))')idx,path(:,idx)
1798 : call wrtout(ount,msg)
1799 : end if
1800 : end do
1801 : end do
1802 320 : path(:,npts)=bounds(:,nbounds)
1803 :
1804 : if (prtvol > 0) then
1805 : write(msg,'(i0,4x,3(f8.5,1x))')npts,path(:,npts)
1806 : call wrtout(ount,msg)
1807 : end if
1808 :
1809 80 : end subroutine make_path
1810 : !!***
1811 :
1812 : !----------------------------------------------------------------------
1813 :
1814 : !!****f* m_bz_mesh/bz_mesh_find_qmesh
1815 : !! NAME
1816 : !! bz_mesh_find_qmesh
1817 : !!
1818 : !! FUNCTION
1819 : !! Find the q-mesh defined as all the possible differences between k-points
1820 : !! Find the irreducible q-points using a special treatment for the Gamma point.
1821 : !! Then call setup_kmesh to initialize the Qmesh datatype.
1822 : !!
1823 : !! INPUTS
1824 : !! Cryst<crystal_t>=datatype gathering info on the unit cell and symmetries
1825 : !! Kmesh<kmesh_t>=datatype gathering information on the k-mesh
1826 : !!
1827 : !! OUTPUT
1828 : !! Qmesh<kmesh_t>=datatype gathering information on the q-point sampling.
1829 : !!
1830 : !! SOURCE
1831 :
1832 221 : subroutine bz_mesh_find_qmesh(Qmesh, Cryst, Kmesh)
1833 :
1834 : !Arguments ------------------------------------
1835 : class(kmesh_t),intent(inout) :: Qmesh
1836 : type(kmesh_t),intent(in) :: Kmesh
1837 : type(crystal_t),intent(in) :: Cryst
1838 :
1839 : !Local variables-------------------------------
1840 : !scalars
1841 : integer :: nqibz, kptopt
1842 : !arrays
1843 221 : real(dp),allocatable :: qibz(:,:)
1844 : ! *************************************************************************
1845 :
1846 : ! Find the number of q-points such that q = k1 - k2.
1847 221 : call findnq(Kmesh%nbz, Kmesh%bz, Cryst%nsym, Cryst%symrec, Cryst%symafm, nqibz, Cryst%timrev)
1848 :
1849 : ! Find the coordinates of the q-points in the IBZ.
1850 663 : ABI_MALLOC(qibz, (3, nqibz))
1851 221 : call findq(Kmesh%nbz, Kmesh%bz, Cryst%nsym, Cryst%symrec, Cryst%symafm, Cryst%gprimd, nqibz, qibz, Cryst%timrev)
1852 :
1853 : ! Create the qmesh object starting from the IBZ.
1854 221 : kptopt = Kmesh%kptopt
1855 221 : call qmesh%init(cryst, nqibz, qibz, kptopt)
1856 221 : ABI_FREE(qibz)
1857 :
1858 221 : end subroutine bz_mesh_find_qmesh
1859 : !!***
1860 :
1861 : !----------------------------------------------------------------------
1862 :
1863 : !!****f* m_bz_mesh/findnq
1864 : !! NAME
1865 : !! findnq
1866 : !!
1867 : !! FUNCTION
1868 : !! Identify the number of q-points in the IBZ by which the k-points in BZ differ
1869 : !! (count the q points in the k-point difference set)
1870 : !!
1871 : !! INPUTS
1872 : !! nkbz=number of k points in Brillouin zone
1873 : !! kbz(3,nkbz)=coordinates of k points in BZ
1874 : !! nsym=number of symmetry operations
1875 : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
1876 : !! symafm(nsym)=-1 if AFM symmetry.
1877 : !! timrev=2 if time-reversal symmetry is used, 1 otherwise
1878 : !!
1879 : !! OUTPUT
1880 : !! nqibz=number of q points
1881 : !!
1882 : !! SOURCE
1883 :
1884 221 : subroutine findnq(nkbz, kbz, nsym, symrec, symafm, nqibz, timrev)
1885 :
1886 : !Arguments ------------------------------------
1887 : !scalars
1888 : integer,intent(in) :: timrev,nkbz,nsym
1889 : integer,intent(out) :: nqibz
1890 : !arrays
1891 : integer,intent(in) :: symrec(3,3,nsym),symafm(nsym)
1892 : real(dp),intent(in) :: kbz(3,nkbz)
1893 :
1894 : !Local variables ------------------------------
1895 : !scalars
1896 : integer :: ifound,ik,isym,iq,memory_exhausted,nqall,nqallm,itim,ierr
1897 : !arrays
1898 : integer :: g0(3)
1899 : real(dp) :: qposs(3),qrot(3)
1900 221 : real(dp),allocatable :: qall(:,:)
1901 : !************************************************************************
1902 :
1903 : ! Infinite do-loop to be able to allocate sufficient memory
1904 221 : nqallm=1000
1905 : do
1906 221 : memory_exhausted=0
1907 663 : ABI_MALLOC_OR_DIE(qall,(3,nqallm), ierr)
1908 : nqall=0
1909 :
1910 : ! Loop over all k-points in BZ, forming k-k1.
1911 10399 : do ik=1,nkbz
1912 40712 : qposs(:)=kbz(:,ik)-kbz(:,1)
1913 :
1914 : ! Check whether this q (or its equivalent) has already been found within a reciprocal lattice vector.
1915 : ! Use spatial inversion instead of time reversal whenever possible.
1916 : ifound=0
1917 90886 : do iq=1,nqall
1918 251966 : do itim=1,timrev
1919 7855972 : do isym=1,nsym
1920 7614184 : if (symafm(isym)==-1) CYCLE
1921 209455456 : qrot = (3-2*itim) * MATMUL(symrec(:,:,isym),qall(:,iq))
1922 7641632 : if (isamek(qrot,qposs,g0)) ifound=ifound+1
1923 : end do
1924 : end do
1925 : end do
1926 :
1927 10399 : if (ifound==0) then
1928 1383 : nqall=nqall+1
1929 :
1930 : ! If not yet found, check that the allocation is big enough.
1931 1383 : if (nqall>nqallm) then
1932 0 : memory_exhausted=1
1933 0 : ABI_FREE(qall)
1934 0 : nqallm=nqallm*2; EXIT ! Exit the do ik=1 loop
1935 : end if
1936 :
1937 : ! Add it to the list.
1938 5532 : qall(:,nqall)=qposs(:)
1939 : end if
1940 : end do
1941 :
1942 221 : if (memory_exhausted==0) EXIT
1943 : end do !infinite loop
1944 :
1945 221 : ABI_FREE(qall)
1946 221 : nqibz=nqall
1947 :
1948 221 : end subroutine findnq
1949 : !!***
1950 :
1951 : !----------------------------------------------------------------------
1952 :
1953 : !!****f* m_bz_mesh/findq
1954 : !! NAME
1955 : !! findq
1956 : !!
1957 : !! FUNCTION
1958 : !! Identify the q-points by which the k-points in BZ differ
1959 : !!
1960 : !! INPUTS
1961 : !! nkbz=number of k points in Brillouin zone
1962 : !! kbz(3,nkbz)=coordinates of k points in BZ
1963 : !! nsym=number of symmetry operations
1964 : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
1965 : !! symafm(nsym)=-1 is symmetry is AFM, +1 otherwise.
1966 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1967 : !! nqibz=number of q points in the IBZ by which k points differ (computed in findnq)
1968 : !! timrev=2 if time-reversal symmetry is used, 1 otherwise
1969 : !!
1970 : !! OUTPUT
1971 : !! qibz(3,nqibz)=coordinates of q points by which k points differ
1972 : !!
1973 : !! SOURCE
1974 :
1975 221 : subroutine findq(nkbz, kbz, nsym, symrec, symafm, gprimd, nqibz, qibz, timrev)
1976 :
1977 : !Arguments ------------------------------------
1978 : !scalars
1979 : integer,intent(in) :: nkbz,nqibz,nsym,timrev
1980 : !arrays
1981 : integer,intent(in) :: symrec(3,3,nsym),symafm(nsym)
1982 : real(dp),intent(in) :: gprimd(3,3),kbz(3,nkbz)
1983 : real(dp),intent(inout) :: qibz(3,nqibz)
1984 :
1985 : !Local variables ------------------------------
1986 : !scalars
1987 : integer :: ii,ik,iq,iqp,isym,itim
1988 : real(dp) :: shift1,qred
1989 : logical :: found
1990 : character(len=500) :: msg
1991 : !arrays
1992 : integer :: g0(3)
1993 : real(dp) :: gmet(3,3),qposs(3),qrot(3)
1994 : !************************************************************************
1995 :
1996 : ! Compute reciprocal space metrics
1997 221 : do ii=1,3
1998 : gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+&
1999 : gprimd(2,ii)*gprimd(2,:)+&
2000 : gprimd(3,ii)*gprimd(3,:)
2001 : end do
2002 : !
2003 : ! Loop over k-points in BZ, form k-k1 and translate in first BZ.
2004 : ! iq is the no. of q-points found, zero at the beginning
2005 221 : iq=0
2006 10399 : do ik=1,nkbz
2007 40712 : qposs(:)=kbz(:,ik)-kbz(:,1)
2008 : ! Check whether this q (or its equivalent) has already been found.
2009 : ! Use spatial inversion instead of time reversal whenever possible.
2010 10178 : found=.FALSE.
2011 90886 : do iqp=1,iq
2012 251966 : do itim=1,timrev
2013 7855972 : do isym=1,nsym
2014 7614184 : if (symafm(isym)==-1) CYCLE
2015 209455456 : qrot = (3-2*itim) * MATMUL(symrec(:,:,isym),qibz(:,iqp))
2016 7641632 : if (isamek(qrot,qposs,g0)) found=.TRUE.
2017 : end do
2018 : end do
2019 : end do
2020 10399 : if (.not.found) then
2021 1383 : iq=iq+1
2022 1383 : if (iq>nqibz) then
2023 0 : ABI_BUG(sjoin('iq > nqibz= ',itoa(nqibz)))
2024 : end if
2025 5532 : qibz(:,iq)=qposs(:)
2026 : end if
2027 : end do
2028 :
2029 221 : if (iq/=nqibz) then
2030 0 : write(msg,'(2(a,i0))')' iq= ',iq,'/= nqibz= ',nqibz
2031 0 : ABI_BUG(msg)
2032 : end if
2033 : !
2034 : ! Translate q-points to 1st BZ in the interval [-1/2,1/2[
2035 1604 : do iq=1,nqibz
2036 5753 : do ii=1,3
2037 4149 : call wrap2_pmhalf(qibz(ii,iq),qred,shift1)
2038 5532 : qibz(ii,iq)=qred
2039 : end do
2040 : end do
2041 :
2042 221 : end subroutine findq
2043 : !!***
2044 :
2045 : !----------------------------------------------------------------------
2046 :
2047 : !!****f* m_bz_mesh/findqg0
2048 : !! NAME
2049 : !! findqg0
2050 : !!
2051 : !! FUNCTION
2052 : !! Identify q + g0 = k - kp
2053 : !!
2054 : !! INPUTS
2055 : !! kmkp(3)= k - kp input vector
2056 : !! nqbz=number of q points in the BZ
2057 : !! qbz(3,nqbz)=coordinates of q-points in the BZ
2058 : !! mG0(3)= For each reduced direction gives the maximum G0 component to account for umklapp processes
2059 : !!
2060 : !! OUTPUT
2061 : !! iq=index of q in qbz array.
2062 : !! g0(3)=reciprocal space vector, to be used in igfft
2063 : !!
2064 : !! SOURCE
2065 :
2066 85293 : subroutine findqg0(iq, g0, kmkp, nqbz, qbz, mG0)
2067 :
2068 : !Arguments ------------------------------------
2069 : !scalars
2070 : integer,intent(in) :: nqbz
2071 : integer,intent(out) :: iq
2072 : !arrays
2073 : integer,intent(in) :: mG0(3)
2074 : integer,intent(out) :: g0(3)
2075 : real(dp),intent(in) :: kmkp(3),qbz(3,nqbz)
2076 :
2077 : !Local variables-------------------------------
2078 : !scalars
2079 : integer :: ig,iqbz,jg01,jg02,jg03
2080 : real(dp) :: tolq0=1.0D-3 !FIXME if I use 1.0d-4 the jobs crash, should understand why
2081 : !character(len=500) :: msg
2082 : !arrays
2083 83104 : real(dp) :: glist1(2*ABS(mG0(1))+1),glist2(2*ABS(mG0(2))+1),glist3(2*ABS(mG0(3))+1), qpg0(3),rg(3)
2084 : ! *************************************************************************
2085 :
2086 85293 : iq = 0
2087 :
2088 110052 : if (ALL(ABS(kmkp) < EPSILON(one))) then
2089 : ! Find q close to 0
2090 121159 : do iqbz=1,nqbz
2091 158434 : if (ALL(ABS(qbz(:,iqbz)) < tolq0)) then
2092 2189 : iq = iqbz
2093 : end if
2094 : end do
2095 :
2096 2189 : ABI_CHECK(iq /= 0, 'Wrong list of q-points: q=0 not present.')
2097 2189 : g0(:) = 0; RETURN
2098 :
2099 : else
2100 : ! q is not zero, find q such as k-kp=q+G0.
2101 :
2102 : ! Try with G0 = 0 first.
2103 : !do iqbz=1,nqbz
2104 : ! if (ALL(ABS(qbz(:,iqbz)-kmkp)<TOL_KDIFF)) then
2105 : ! iq=iqbz
2106 : ! g0(:)=0; RETURN
2107 : ! end if
2108 : !end do
2109 :
2110 : ! Init G0 lists to accelerate search below (small |G0| first)
2111 83104 : glist1(1) = 0; ig = 2
2112 169424 : do jg01=1,mG0(1)
2113 86320 : glist1(ig) = jg01
2114 86320 : glist1(ig+1) = -jg01
2115 169424 : ig = ig + 2
2116 : end do
2117 :
2118 83104 : glist2(1) = 0; ig = 2
2119 210152 : do jg02=1,mG0(2)
2120 127048 : glist2(ig) = jg02
2121 127048 : glist2(ig+1) = -jg02
2122 210152 : ig = ig + 2
2123 : end do
2124 :
2125 83104 : glist3(1) = 0; ig = 2
2126 169668 : do jg03=1,mG0(3)
2127 86564 : glist3(ig) = jg03
2128 86564 : glist3(ig+1) = -jg03
2129 169668 : ig = ig + 2
2130 : end do
2131 :
2132 : ! Search algorithm.
2133 113965 : g1loop: do jg01=1,2*mG0(1)+1
2134 113965 : rg(1) = glist1(jg01)
2135 272590 : do jg02=1,2*mG0(2)+1
2136 241729 : rg(2) = glist2(jg02)
2137 794783 : do jg03=1,2*mG0(3)+1
2138 605297 : rg(3) = glist3(jg03)
2139 :
2140 : ! Form q+G0 and check if it is the one.
2141 38796478 : do iqbz=1,nqbz
2142 152462640 : qpg0= qbz(:,iqbz) + rg
2143 42060906 : if (ALL(ABS(qpg0-kmkp) < TOL_KDIFF)) then
2144 83104 : iq = iqbz
2145 332416 : g0 = NINT(rg)
2146 : EXIT g1loop
2147 : end if
2148 : end do
2149 :
2150 : end do
2151 : end do
2152 : end do g1loop
2153 :
2154 83104 : if (iq == 0) then
2155 0 : ABI_ERROR(sjoin('q = k-kp+G0 not found. kmkp:', ktoa(kmkp)))
2156 : end if
2157 : end if
2158 :
2159 : end subroutine findqg0
2160 : !!***
2161 :
2162 : !----------------------------------------------------------------------
2163 :
2164 : !!****f* m_bz_mesh/littlegroup_init
2165 : !! NAME
2166 : !! littlegroup_init
2167 : !!
2168 : !! FUNCTION
2169 : !! Finds symmetry operations belonging to the little group associated to an external
2170 : !! point ext_pt and fills symmetry tables.
2171 : !!
2172 : !! INPUTS
2173 : !! ext_pt(3)= External point in the Brillouin zone in reduce coordinated
2174 : !! nbz=number of points in the full BZ.
2175 : !! bz(3,nbz)=points in the full BZ.
2176 : !! Cryst<crystal_t>= Info on symmetries and unit cell.
2177 : !! use_umklp=flag to include umklapp G0 vectors in the definition of the little group (0:n0,1:yes)
2178 : !! npwe=If greater than 0, the index of G-Go in the gvec(:,1:npwvec) array will be calculated
2179 : !! and stored in %igmG0 for each symmetry preserving the external q. Note that G is one of the npwe vectors.
2180 : !! gvec(3,npwe) coordinates of G vectors
2181 : !! [timrev]=Optional argument to change the value of time-reversal. If not given, the value from cryst is used.
2182 : !!
2183 : !! SOURCE
2184 :
2185 1052 : subroutine littlegroup_init(Ltg, ext_pt, nbz, bz, Cryst, use_umklp, npwe, gvec, timrev)
2186 :
2187 : !Arguments ------------------------------------
2188 : !scalars
2189 : class(littlegroup_t),intent(inout) :: Ltg
2190 : integer,intent(in) :: nbz, npwe, use_umklp
2191 : type(crystal_t),target,intent(in) :: Cryst
2192 : real(dp),intent(in) :: bz(3, nbz)
2193 : !arrays
2194 : integer,optional,intent(in) :: gvec(:,:) ! (3,npwe)
2195 : real(dp),intent(in) :: ext_pt(3)
2196 : integer,optional,intent(in) :: timrev
2197 :
2198 : !Local variables-------------------------------
2199 : !scalars
2200 : integer :: dummy_timrev,enough,idx,ige,igpw,ik,ind,iold,iout,isym,itest,itim
2201 : integer :: nkibzq,nsym,nsym_Ltg,ntest,my_timrev,ierr,npwvec
2202 : real(dp) :: G0len,kin,mG0len,max_kin
2203 : logical :: found,found_identity,use_antiferro
2204 : character(len=500) :: msg
2205 : !arrays
2206 2104 : integer :: g0(3),gg(3),gmG0(3),identity(3,3),nop(Cryst%timrev),nopg0(2)
2207 2104 : integer :: symxpt(4,2,Cryst%nsym)
2208 1052 : integer,allocatable :: indkpt1(:),symafm_ltg(:),symrec_Ltg(:,:,:),bz2ibz_smap(:,:)
2209 1052 : integer,pointer :: symafm(:),symrec(:,:,:)
2210 : real(dp) :: knew(3)
2211 1052 : real(dp),allocatable :: ktest(:,:),wtk(:),wtk_folded(:)
2212 : !************************************************************************
2213 :
2214 2098 : ABI_CHECK(any(cryst%timrev == [1, 2]), sjoin("Wrong value for cryst%timrev:", itoa(cryst%timrev)))
2215 :
2216 : ! Destroy structure if it already exists
2217 1052 : call Ltg%free()
2218 :
2219 : ! Copy useful data.
2220 1052 : nsym = Cryst%nsym
2221 1052 : my_timrev = Cryst%timrev
2222 1052 : if (present(timrev)) then
2223 0 : my_timrev = timrev
2224 0 : ABI_CHECK_ILEQ(my_timrev, cryst%timrev, "my_timrev cannot be greater that cryst%timrev")
2225 : end if
2226 :
2227 1052 : symrec => Cryst%symrec
2228 1052 : symafm => Cryst%symafm
2229 1052 : use_antiferro = Cryst%use_antiferro
2230 :
2231 : ! Store dimensions and useful info.
2232 1052 : Ltg%nsym_sg = nsym
2233 1052 : Ltg%timrev = my_timrev
2234 1052 : Ltg%nbz = nbz
2235 : !Ltg%use_umklp=use_umklp ! 1 if umklapp processes are used
2236 4208 : Ltg%ext_pt(:)=ext_pt(:)
2237 :
2238 3156 : ABI_MALLOC(Ltg%G0, (3,2,nsym))
2239 3156 : ABI_MALLOC(Ltg%ibzq, (nbz))
2240 2104 : ABI_MALLOC(Ltg%bz2ibz, (nbz))
2241 3156 : ABI_MALLOC(Ltg%preserve, (2, nsym))
2242 4208 : ABI_MALLOC(Ltg%wtksym, (2, nsym, nbz))
2243 2104 : ABI_MALLOC(Ltg%tab, (nbz))
2244 2104 : ABI_MALLOC(Ltg%tabi, (nbz))
2245 2104 : ABI_MALLOC(Ltg%tabo, (nbz))
2246 2104 : ABI_MALLOC(Ltg%flag_umklp, (2, nsym))
2247 :
2248 : ! In the old GW implementation we were removing symmetries related by time-reversal and
2249 : ! sometimes it happened that only the inversion was reported in the KSS file (see outkss.F90).
2250 1052 : identity(:,:)=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/)) ; found_identity=.FALSE.
2251 1052 : do isym=1,nsym
2252 13676 : if (ALL(symrec(:,:,isym) == identity)) then
2253 : found_identity=.TRUE.; EXIT
2254 : end if
2255 : end do
2256 1052 : if (.not. found_identity) then
2257 : write(msg,'(5a)')&
2258 0 : 'Only the inversion was found in the set of symmetries read from the KSS file ',ch10,&
2259 0 : 'Likely you are using a KSS file generated with an old version of Abinit, ',ch10,&
2260 0 : 'To run a GW calculation with an old KSS file, use version < 5.5 '
2261 0 : ABI_ERROR(msg)
2262 : end if
2263 :
2264 : ! Find operations in the little group as well as umklapp vectors G0
2265 1052 : call littlegroup_q(nsym, ext_pt, symxpt, symrec, symafm, dummy_timrev, prtvol=0)
2266 :
2267 703412 : Ltg%preserve(:,:)=0; Ltg%g0(:,:,:)=0; Ltg%flag_umklp(:,:)=0; mG0len=zero
2268 :
2269 3150 : do itim=1,my_timrev
2270 96510 : do isym=1,nsym
2271 93360 : if (symafm(isym)==-1) CYCLE
2272 :
2273 92002 : if (symxpt(4, itim, isym) == 1) then !\pm Sq = q+g0
2274 108227 : if (ANY(symxpt(1:3, itim, isym) /= 0) .and. use_umklp == 0) CYCLE ! Exclude non zero G0 vectors
2275 25781 : Ltg%preserve(itim, isym) = 1
2276 180467 : g0(:)=symxpt(1:3, itim, isym); Ltg%g0(:, itim, isym) = g0(:)
2277 101930 : if (ANY(Ltg%g0(:, itim, isym) /= 0)) Ltg%flag_umklp(itim, isym) = 1
2278 : ! Max radius to be considered to include all G0s
2279 25781 : G0len = normv(g0,Cryst%gmet,'G')
2280 : mG0len = MAX(mG0len, G0len)
2281 : end if
2282 : end do
2283 : end do
2284 :
2285 3150 : nop(:) = 0; nopg0(:) = 0
2286 3150 : do itim=1,my_timrev
2287 95458 : nop (itim) = SUM(Ltg%preserve (itim,:))
2288 1052 : nopg0(itim) = SUM(Ltg%flag_umklp(itim,:))
2289 : end do
2290 3150 : nsym_Ltg = SUM(nop(:))
2291 :
2292 : ! Store little group operations, include time-reversal if present.
2293 1052 : Ltg%nsym_Ltg = nsym_Ltg
2294 3156 : ABI_MALLOC(symrec_Ltg, (3, 3, Ltg%nsym_Ltg))
2295 :
2296 1052 : ind = 1
2297 3150 : do itim=1,my_timrev
2298 96510 : do isym=1,nsym
2299 95458 : if (Ltg%preserve(itim,isym)==1) then
2300 182225 : if (itim==1) symrec_Ltg(:,:,ind) = symrec(:,:,isym)
2301 178709 : if (itim==2) symrec_Ltg(:,:,ind) =-symrec(:,:,isym)
2302 25781 : ind = ind+1
2303 : end if
2304 : end do
2305 : end do
2306 :
2307 : ! Check the closure of the (ferromagnetic) little group
2308 3156 : ABI_MALLOC(symafm_ltg,(Ltg%nsym_Ltg))
2309 26833 : symafm_ltg(:) = 1
2310 1052 : call sg_multable(Ltg%nsym_Ltg,symafm_ltg,symrec_Ltg,ierr)
2311 1052 : ABI_CHECK(ierr == 0, "Error in group closure")
2312 :
2313 1052 : ABI_FREE(symafm_ltg)
2314 :
2315 : ! Find the irreducible zone associated to ext_pt
2316 : ! Do not use time-reversal since it has been manually introduced previously
2317 2104 : ABI_MALLOC(indkpt1, (nbz))
2318 3156 : ABI_MALLOC(wtk_folded, (nbz))
2319 2104 : ABI_MALLOC(wtk, (nbz))
2320 3156 : ABI_MALLOC(bz2ibz_smap, (6, nbz))
2321 85926 : wtk=one; iout=0; dummy_timrev=0
2322 :
2323 : call symkpt(0,Cryst%gmet,indkpt1,iout, bz, nbz, nkibzq, Ltg%nsym_Ltg, symrec_Ltg, dummy_timrev, wtk, wtk_folded, &
2324 1052 : bz2ibz_smap, xmpi_comm_self)
2325 :
2326 1052 : ABI_FREE(bz2ibz_smap)
2327 1052 : ABI_FREE(indkpt1)
2328 1052 : ABI_FREE(wtk)
2329 :
2330 1052 : Ltg%nibz_Ltg = nkibzq
2331 : !
2332 : ! === Set up table in the BZ ===
2333 : ! * 0 if the point does not belong to IBZ_xpt, 1 otherwise
2334 3156 : ABI_MALLOC(Ltg%ibz2bz,(nkibzq))
2335 204021 : Ltg%ibzq(:)=0; Ltg%bz2ibz(:)=0; Ltg%ibz2bz(:)=0
2336 :
2337 1052 : ind=0; enough=0
2338 85926 : do ik=1,nbz
2339 85926 : if (wtk_folded(ik)>tol8) then
2340 33221 : ind = ind + 1
2341 33221 : Ltg%ibzq(ik) = 1
2342 33221 : Ltg%bz2ibz(ik) = ind
2343 33221 : Ltg%ibz2bz(ind)= ik
2344 : end if
2345 : end do
2346 1052 : ABI_CHECK_IEQ(ind, Ltg%nibz_Ltg, "BUG: ind /= Ltg%nibz_Ltg")
2347 :
2348 : ! Reconstruct full BZ starting from IBZ_q.
2349 : ! Calculate appropriate weight for each item (point,symmetry operation,time-reversal)
2350 12376956 : Ltg%tab=0; Ltg%tabo=0; Ltg%tabi=0; Ltg%wtksym(:,:,:)=0
2351 :
2352 : ! Start with zero no. of k-points found
2353 1052 : ntest = 0
2354 3156 : ABI_MALLOC(ktest, (3, nbz))
2355 340548 : ktest = zero
2356 :
2357 85926 : do ik=1,nbz
2358 84874 : if (Ltg%ibzq(ik) /= 1) CYCLE
2359 : ! Loop over symmetry operations S and time-reversal.
2360 : ! Use spatial inversion instead of time reversal whenever possible.
2361 100693 : do itim=1,my_timrev
2362 3281742 : do isym=1,nsym
2363 :
2364 : ! Form IS k only for (IS) pairs in the (ferromagnetic) little group.
2365 3130448 : if (symafm(isym)==-1) CYCLE
2366 3051728 : if (Ltg%preserve(itim,isym)==0) CYCLE
2367 9651600 : knew(:)=(3-2*itim)*MATMUL(symrec(:,:,isym), bz(:,ik))
2368 : !
2369 : ! Check whether it has already been found (to within a RL vector)
2370 : iold=0
2371 26476705 : do itest=1,ntest
2372 26476705 : if (isamek(knew(:),ktest(:,itest),gg)) iold=iold+1
2373 : end do
2374 :
2375 411120 : if (iold==0) then
2376 : ! Found new BZ point
2377 : ! For this point the operation (isym,itim) must be considered to reconstruct the full BZ
2378 84874 : Ltg%wtksym(itim,isym,ik)=1
2379 84874 : ntest=ntest+1
2380 339496 : ktest(:,ntest)=knew(:)
2381 : !
2382 : ! Now find knew in the BZ array
2383 7258112 : found=.FALSE.
2384 7258112 : do idx=1,nbz
2385 7258112 : if (isamek(knew(:), bz(:,idx), gg)) then ! They are the same within a RL vector
2386 84874 : Ltg%tab (idx)=ik
2387 84874 : Ltg%tabo(idx)=isym
2388 84874 : Ltg%tabi(idx)=3-2*itim
2389 84874 : found=.TRUE.; EXIT
2390 : end if
2391 : end do
2392 : if (.not.found) then
2393 0 : write(msg,'(a,3f12.6,a)')'Not able to find the ',knew(:),' in the array BZ '
2394 0 : ABI_ERROR(msg)
2395 : end if
2396 : end if
2397 :
2398 : end do ! isym
2399 : end do ! itim
2400 : end do ! nbz
2401 :
2402 1052 : ABI_FREE(ktest)
2403 :
2404 1052 : if (ntest/=nbz) then
2405 0 : ABI_BUG(sjoin('ntest - nbz = ',itoa(ntest-nbz)))
2406 : end if
2407 :
2408 12122334 : if (sum(Ltg%wtksym) /= nbz) then
2409 0 : ABI_BUG(sjoin('sum(Ltg%wtksym)-nbz = ', itoa(SUM(Ltg%wtksym)-nbz)))
2410 : end if
2411 :
2412 1052 : Ltg%max_kin_gmG0=zero
2413 :
2414 1052 : if (npwe > 0.and. PRESENT(gvec)) then
2415 934 : npwvec = SIZE(gvec,DIM=2)
2416 : ! This correspond to the case in which we need to know the index of G-Go in the gvec array
2417 : ! where G is one of the npwe vectors. This is required in screening but not in sigma.
2418 : ! The drawback is that the effective G sphere used to calculate the oscillators must be smaller
2419 : ! that gvec if we want to avoid possible aliasing effects. Lifting this constraint would require
2420 : ! a lot of boring coding. (no need to do this if ext_pt=zero, but oh well)
2421 3736 : ABI_MALLOC(Ltg%igmG0,(npwe,2,nsym))
2422 4604494 : Ltg%igmG0(:,:,:)=0
2423 : max_kin=zero
2424 :
2425 : ! Loop over symmetry operations S and time-reversal
2426 2796 : do itim=1,my_timrev
2427 85468 : do isym=1,nsym
2428 : ! Form IS k only for (IS) pairs in the little group
2429 82672 : if (symafm(isym)==-1) CYCLE
2430 81654 : if (Ltg%preserve(itim,isym)/=0) then
2431 76564 : g0(:)=Ltg%g0(:,itim,isym)
2432 1107614 : do ige=1,npwe
2433 4353892 : gmG0(:)=gvec(:,ige)-g0(:)
2434 1088473 : kin=half*normv(gmG0,Cryst%gmet,'G')**2
2435 1088473 : max_kin=MAX(max_kin,kin)
2436 :
2437 1088473 : found=.FALSE.
2438 55354615 : do igpw=1,npwvec
2439 69332209 : if (ALL(gvec(:,igpw)-gmG0(:)==0)) then
2440 1088473 : Ltg%igmG0(ige,itim,isym)=igpw
2441 : found=.TRUE.; EXIT
2442 : end if
2443 : end do
2444 19141 : if (.not. found) then
2445 : write(msg,'(5a,f8.3,2a,3i5)')&
2446 0 : 'Not able to found G-G0 in the largest G-spere ',ch10,&
2447 0 : 'Decrease the size of epsilon or, if possible, increase ecutwfn (>ecuteps) ',ch10,&
2448 0 : 'Minimum required cutoff energy for G-G0 sphere= ',kin,ch10,&
2449 0 : 'G0 = ',g0(:)
2450 0 : ABI_ERROR(msg)
2451 : end if
2452 : end do
2453 : end if
2454 : end do
2455 : end do
2456 934 : Ltg%max_kin_gmG0=max_kin
2457 : end if
2458 1052 : ABI_FREE(symrec_Ltg)
2459 :
2460 : #if 0
2461 : ! DEBUGGING SECTION
2462 : do ik=1,nbz
2463 : if (ABS(SUM(Ltg%wtksym(1,:,ik)+Ltg%wtksym(2,:,ik))-wtk_folded(ik))>tol6) then
2464 : write(std_out,*)' sum(Ltg%wtksym,ik)-wtk_folded(ik) = ',sum(Ltg%wtksym(1,:,ik)+Ltg%wtksym(2,:,ik))-wtk_folded(ik)
2465 : write(std_out,*)Ltg%wtksym(1,:,ik),Ltg%wtksym(2,:,ik),wtk_folded(ik)
2466 : write(std_out,*)ik, bz(:,ik)
2467 : ABI_BUG("Wrong weight")
2468 : end if
2469 : end do
2470 : do ik=1,nbz
2471 : knew = Ltg%tabi(ik) * MATMUL(symrec(:,:,Ltg%tabo(ik)), bz(:,Ltg%tab(ik)))
2472 : if (.not.isamek(knew, bz(:,ik),gg)) then
2473 : write(std_out,*)knew, bz(:,ik)
2474 : write(std_out,*)Ltg%tabo(ik),Ltg%tabi(ik),Ltg%tab(ik)
2475 : ABI_BUG("Wrong tables")
2476 : end if
2477 : end do
2478 : #endif
2479 :
2480 1052 : ABI_FREE(wtk_folded)
2481 :
2482 : DBG_EXIT("COLL")
2483 :
2484 2104 : end subroutine littlegroup_init
2485 : !!***
2486 :
2487 : !----------------------------------------------------------------------
2488 :
2489 : !!****f* m_bz_mesh/littlegroup_free_0D
2490 : !! NAME
2491 : !! littlegroup_free_0D
2492 : !!
2493 : !! FUNCTION
2494 : !! Deallocate dynamic memory
2495 : !!
2496 : !! SOURCE
2497 :
2498 2625 : subroutine littlegroup_free_0D(Ltg)
2499 :
2500 : !Arguments ------------------------------------
2501 : class(littlegroup_t),intent(inout) :: Ltg
2502 : ! *********************************************************************
2503 :
2504 2625 : ABI_SFREE(Ltg%g0)
2505 2625 : ABI_SFREE(Ltg%ibzq)
2506 2625 : ABI_SFREE(Ltg%bz2ibz)
2507 2625 : ABI_SFREE(Ltg%ibz2bz)
2508 2625 : ABI_SFREE(Ltg%igmG0)
2509 2625 : ABI_SFREE(Ltg%flag_umklp)
2510 2625 : ABI_SFREE(Ltg%preserve)
2511 2625 : ABI_SFREE(Ltg%tab)
2512 2625 : ABI_SFREE(Ltg%tabo)
2513 2625 : ABI_SFREE(Ltg%tabi)
2514 2625 : ABI_SFREE(Ltg%wtksym)
2515 :
2516 2625 : end subroutine littlegroup_free_0D
2517 : !!***
2518 :
2519 : !----------------------------------------------------------------------
2520 :
2521 : !!****f* m_bz_mesh/littlegroup_free_1D
2522 : !! NAME
2523 : !! littlegroup_free_1D
2524 : !!
2525 : !! FUNCTION
2526 : !! Deallocate dynamic memory
2527 : !!
2528 : !! SOURCE
2529 :
2530 349 : subroutine littlegroup_free_1D(Ltg)
2531 :
2532 : !Arguments ------------------------------------
2533 : class(littlegroup_t),intent(inout) :: Ltg(:)
2534 :
2535 : !Local variables-------------------------------
2536 : integer :: ipt
2537 : ! *********************************************************************
2538 :
2539 1922 : do ipt=1,SIZE(Ltg)
2540 1922 : call littlegroup_free_0D(Ltg(ipt))
2541 : end do
2542 :
2543 349 : end subroutine littlegroup_free_1D
2544 : !!***
2545 :
2546 : !----------------------------------------------------------------------
2547 :
2548 : !!****f* m_bz_mesh/littlegroup_print
2549 : !! NAME
2550 : !! littlegroup_print
2551 : !!
2552 : !! FUNCTION
2553 : !! Print info on the littlegroup_t data type.
2554 : !!
2555 : !! INPUTS
2556 : !! Ltg=the datatype to be printed
2557 : !! units=unit numbers for output
2558 : !! [prtvol]=verbosity level
2559 : !!
2560 : !! OUTPUT
2561 : !! Only printing
2562 : !!
2563 : !! SOURCE
2564 :
2565 1061 : subroutine littlegroup_print(Ltg, units, prtvol)
2566 :
2567 : !Arguments ------------------------------------
2568 : class(littlegroup_t),intent(in) :: Ltg
2569 : integer,intent(in) :: units(:)
2570 : integer,optional,intent(in) :: prtvol
2571 :
2572 : !Local variables-------------------------------
2573 : !scalars
2574 : integer :: itim, my_prtvol
2575 : character(len=500) :: msg
2576 : !arrays
2577 2122 : integer :: nop(Ltg%timrev),nopg0(Ltg%timrev)
2578 : ! *********************************************************************
2579 :
2580 1061 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
2581 :
2582 1061 : write(msg,'(7a,i0,a,i0,2a,i0,a,i0)')ch10, &
2583 1061 : ' ==== Little Group Info ==== ', ch10, &
2584 1061 : ' External point: ',trim(ktoa(Ltg%ext_pt)), ch10, &
2585 1061 : ' Number of points in the IBZ defined by little group: ', Ltg%nibz_Ltg, '/', Ltg%nbz,ch10, &
2586 2122 : ' Number of operations in the little group: ',Ltg%nsym_Ltg,'/',Ltg%nsym_sg
2587 1061 : call wrtout(units, msg)
2588 :
2589 6360 : nop=0 ; nopg0=0
2590 3180 : do itim=1,Ltg%timrev
2591 96679 : nop (itim)=SUM(Ltg%preserve (itim,:))
2592 97740 : nopg0(itim)=SUM(Ltg%flag_umklp(itim,:))
2593 : end do
2594 :
2595 3180 : do itim=1,Ltg%timrev
2596 3180 : if (itim==1) then
2597 : write(msg,'(2(a,i2,a))') &
2598 1061 : ' No time-reversal symmetry with zero umklapp: ',nop(1)-nopg0(1),ch10,&
2599 2122 : ' No time-reversal symmetry with non-zero umklapp: ',nopg0(1),ch10
2600 1061 : call wrtout(units, msg)
2601 1058 : else if (itim==2) then
2602 : write(msg,'(2(a,i2,a))') &
2603 1058 : ' time-reversal symmetry with zero umklapp: ',nop(2)-nopg0(2),ch10,&
2604 2116 : ' time-reversal symmetry with non-zero umklapp: ',nopg0(2),ch10
2605 1058 : call wrtout(units, msg)
2606 : end if
2607 : end do
2608 :
2609 1061 : end subroutine littlegroup_print
2610 : !!***
2611 :
2612 : !----------------------------------------------------------------------
2613 :
2614 : !!****f* m_bz_mesh/box_len
2615 : !! NAME
2616 : !! box_len
2617 : !!
2618 : !! FUNCTION
2619 : !! Given a direction in q-space defined by the q-point qpt, this function returns
2620 : !! the length of the vector connecting the origin with one the faces of the cell
2621 : !! defined by the lattice vectors gprimd.
2622 : !!
2623 : !! INPUTS
2624 : !! qpt(3)=The reduced coordinates of the q-point defining the direction. Normalization is not mandatory.
2625 : !! gprimd(3,3)=Cartesian coordinates of the vectors defining the lattice.
2626 : !!
2627 : !! SOURCE
2628 :
2629 0 : function box_len(qpt, gprimd)
2630 :
2631 : !Arguments ------------------------------------
2632 : !scalars
2633 : real(dp) :: box_len
2634 : !arrays
2635 : real(dp),intent(in) :: qpt(3),gprimd(3,3)
2636 :
2637 : !Local variables-------------------------------
2638 : !scalars
2639 : integer :: idir,iplane
2640 : real(dp) :: x1,x2,x3
2641 : !arrays
2642 : real(dp) :: my_qpt(3),gmet(3,3),q0box(3)
2643 : ! *************************************************************************
2644 :
2645 : ! Compute reciprocal space metric
2646 0 : gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
2647 :
2648 : ! Rotate the input q-point such that it is always in the first octant then normalize it.
2649 : ! Bravais lattices are invariant under inversion of any of the basis vectors.
2650 0 : my_qpt = ABS(qpt)/normv(qpt,gmet,"G")
2651 :
2652 : ! Check whether the q is along one of the reduced directions.
2653 0 : idir=0; if (COUNT(ABS(qpt)<tol16) == 2) idir = imax_loc(ABS(qpt))
2654 :
2655 0 : if (idir/=0) then ! easy as q is along vector idir.
2656 0 : box_len = normv(gprimd(:,idir), gmet, "G")
2657 0 : RETURN
2658 :
2659 : else
2660 : !iplane/=0 means that q is placed on one the planes defined by two reciprocal lattice vectors.
2661 0 : iplane=0; if (COUNT(ABS(qpt)<tol16) == 1) iplane = imin_loc(ABS(qpt))
2662 0 : q0box = [-1,-1,-1]
2663 :
2664 0 : if (iplane/=1) then
2665 0 : x1=one
2666 0 : x2=my_qpt(2)/my_qpt(1)
2667 0 : x3=my_qpt(3)/my_qpt(1)
2668 0 : if (x2<=one+tol16 .and. x3<=one+tol16) q0box=(/x1,x2,x3/)
2669 : end if
2670 0 : if (iplane/=2) then
2671 0 : x1=my_qpt(1)/my_qpt(2)
2672 0 : x2=one
2673 0 : x3=my_qpt(3)/my_qpt(2)
2674 0 : if (x1<=one+tol16 .and. x3<=one+tol16) q0box=(/x1,x2,x3/)
2675 : end if
2676 0 : if (iplane/=3) then
2677 0 : x1=my_qpt(1)/my_qpt(3)
2678 0 : x2=my_qpt(2)/my_qpt(3)
2679 0 : x3=one
2680 0 : if (x1<=one+tol16 .and. x2<=one+tol16) q0box=(/x1,x2,x3/)
2681 : end if
2682 :
2683 0 : if (ALL(q0box == [-1,-1,-1])) then
2684 0 : ABI_BUG("Cannot found q0box")
2685 : end if
2686 :
2687 0 : box_len = normv(q0box,gmet,"G")
2688 0 : RETURN
2689 : end if
2690 :
2691 : end function box_len
2692 : !!***
2693 :
2694 : !----------------------------------------------------------------------
2695 :
2696 : !!****f* m_bz_mesh/kpath_init
2697 : !! NAME
2698 : !! kpath_init
2699 : !!
2700 : !! FUNCTION
2701 : !! Create a normalized path given the extrema.
2702 : !!
2703 : !! INPUTS
2704 : !! bounds(3,nbounds)=The points defining the path in reduced coordinates.
2705 : !! gprimd(3,3)=Reciprocal lattice vectors
2706 : !! ndivsm=Number of divisions to be used for the smallest segment.
2707 : !! A negative value activates a specialized mode in which bounds supplies the full list of k-points.
2708 : !!
2709 : !! SOURCE
2710 :
2711 49 : subroutine kpath_init(kpath, bounds, gprimd, ndivsm)
2712 :
2713 : !Arguments ------------------------------------
2714 : !scalars
2715 : class(kpath_t),intent(out) :: kpath
2716 : integer,intent(in) :: ndivsm
2717 : !!arrays
2718 : real(dp),intent(in) :: bounds(:,:),gprimd(3,3)
2719 :
2720 : !Local variables-------------------------------
2721 : integer :: ii
2722 : !arrays
2723 : real(dp) :: dk(3)
2724 : ! *************************************************************************
2725 :
2726 49 : ABI_CHECK(size(bounds, dim=1) == 3, "Wrong dim1 in bounds")
2727 : !ABI_CHECK(ndivsm > 0, sjoin("ndivsm:", itoa(ndivsm)))
2728 49 : Kpath%nbounds = size(bounds, dim=2)
2729 49 : Kpath%ndivsm = ndivsm
2730 :
2731 : ! Compute reciprocal space metric.
2732 3234 : Kpath%gprimd = gprimd; Kpath%gmet = matmul(transpose(gprimd), gprimd)
2733 :
2734 147 : ABI_MALLOC(Kpath%ndivs, (Kpath%nbounds-1))
2735 :
2736 49 : if (kpath%ndivsm > 0) then
2737 47 : call make_path(Kpath%nbounds, bounds, Kpath%gmet, "G", ndivsm, Kpath%ndivs, Kpath%npts, kpath%points, unit=dev_null)
2738 : else
2739 : ! Get list of points directly.
2740 4 : kpath%ndivs = 1
2741 2 : kpath%npts = kpath%nbounds
2742 6 : ABI_MALLOC(Kpath%points, (3, Kpath%npts))
2743 20 : kpath%points = bounds
2744 : end if
2745 :
2746 147 : ABI_MALLOC(Kpath%bounds, (3, Kpath%nbounds))
2747 826 : Kpath%bounds = bounds
2748 :
2749 : ! Compute distance between point i-1 and i
2750 3257 : ABI_CALLOC(kpath%dl, (kpath%npts))
2751 3110 : do ii=2,kpath%npts
2752 12244 : dk = kpath%points(:, ii-1) - kpath%points(:,ii)
2753 3110 : kpath%dl(ii) = normv(dk, kpath%gmet, "G")
2754 : end do
2755 :
2756 : ! Mapping bounds --> points
2757 147 : ABI_MALLOC(kpath%bounds2kpt, (kpath%nbounds))
2758 49 : kpath%bounds2kpt(1) = 1
2759 182 : do ii=1,kpath%nbounds-1
2760 549 : kpath%bounds2kpt(ii+1) = sum(kpath%ndivs(:ii)) + 1
2761 : end do
2762 :
2763 49 : end subroutine kpath_init
2764 : !!***
2765 :
2766 : !----------------------------------------------------------------------
2767 :
2768 : !!****f* m_bz_mesh/kpath_get_versors
2769 : !! NAME
2770 : !! kpath_get_versors
2771 : !!
2772 : !! FUNCTION
2773 : !! Return all the versors emanating from the Gamma point.
2774 : !!
2775 : !! OUTPUT
2776 : !! nvers=number of versors
2777 : !! red_versors(3,nvers)=versors in reduced coords
2778 : !! cart_versors(3,nvers)=versors in reduced coords
2779 : !!
2780 : !! SOURCE
2781 :
2782 1 : subroutine kpath_get_versors(kpath, nvers, red_versors, cart_versors)
2783 :
2784 : !Arguments ------------------------------------
2785 : !scalars
2786 : class(kpath_t),intent(in) :: kpath
2787 : integer,intent(out) :: nvers
2788 : real(dp),allocatable,intent(out) :: red_versors(:,:), cart_versors(:,:)
2789 :
2790 : ! local variables
2791 1 : integer :: ii, ipt, ipt_list(kpath%npts), cnt
2792 : !real(dp) :: norm
2793 1 : real(dp),allocatable :: tmp_versors(:,:)
2794 : ! *************************************************************************
2795 :
2796 : ! Quick return if just one point.
2797 1 : if (kpath%npts == 1) then
2798 0 : nvers = 0
2799 0 : ABI_MALLOC(red_versors, (0, 0))
2800 0 : ABI_MALLOC(cart_versors, (0, 0))
2801 : return
2802 : end if
2803 :
2804 : cnt = 0
2805 43 : do ipt=1,kpath%npts
2806 169 : if (sum(kpath%points(:,ipt)**2) < tol14) then
2807 2 : cnt = cnt + 1
2808 2 : ipt_list(cnt) = ipt
2809 : end if
2810 : end do
2811 :
2812 3 : ABI_MALLOC(tmp_versors, (3, 2*cnt))
2813 1 : nvers = 0
2814 3 : do ii=1,cnt
2815 2 : ipt = ipt_list(ii)
2816 2 : nvers = nvers + 1
2817 : ! Different logic depending whether Gamma is at the beginning/end of the path or in the middle.
2818 3 : if (ipt == 1) then
2819 4 : tmp_versors(:, nvers) = kpath%points(:, ipt+1) - kpath%points(:, ipt)
2820 1 : else if (ipt == kpath%npts) then
2821 4 : tmp_versors(:, nvers) = kpath%points(:, ipt-1) - kpath%points(:, ipt)
2822 : else
2823 0 : tmp_versors(:, nvers) = kpath%points(:, ipt-1) - kpath%points(:, ipt)
2824 0 : nvers = nvers + 1
2825 0 : tmp_versors(:, nvers) = kpath%points(:, ipt+1) - kpath%points(:, ipt)
2826 : end if
2827 : end do
2828 :
2829 : ! Allocate output results
2830 3 : ABI_MALLOC(red_versors, (3, nvers))
2831 10 : red_versors = tmp_versors(:,1:nvers)
2832 1 : ABI_FREE(tmp_versors)
2833 :
2834 : ! Normalize
2835 : !do ii=1,nvers
2836 : ! norm = dot_product(red_versors(:,ii), matmul(kpath%gmet, red_versors(:,ii)))
2837 : ! red_versors(:,ii) = red_versors(:,ii) / sqrt(norm)
2838 : !end do
2839 :
2840 : ! Convert to Cartesian coordinates.
2841 2 : ABI_MALLOC(cart_versors, (3, nvers))
2842 3 : do ii=1, nvers
2843 33 : cart_versors(:,ii) = matmul(kpath%gprimd, red_versors(:,ii))
2844 : end do
2845 :
2846 : end subroutine kpath_get_versors
2847 : !!***
2848 :
2849 : !----------------------------------------------------------------------
2850 :
2851 : !!****f* m_bz_mesh/kpath_free
2852 : !! NAME
2853 : !! kpath_free
2854 : !!
2855 : !! FUNCTION
2856 : !! Free memory allocated in the object
2857 : !!
2858 : !! SOURCE
2859 :
2860 122 : subroutine kpath_free(Kpath)
2861 :
2862 : !Arguments ------------------------------------
2863 : class(kpath_t),intent(inout) :: Kpath
2864 : ! *************************************************************************
2865 :
2866 122 : ABI_SFREE(Kpath%ndivs)
2867 122 : ABI_SFREE(Kpath%bounds2kpt)
2868 122 : ABI_SFREE(Kpath%bounds)
2869 122 : ABI_SFREE(Kpath%points)
2870 122 : ABI_SFREE(Kpath%dl)
2871 :
2872 122 : end subroutine kpath_free
2873 : !!***
2874 :
2875 : !----------------------------------------------------------------------
2876 :
2877 : !!****f* m_bz_mesh/kpath_print
2878 : !! NAME
2879 : !! kpath_print
2880 : !!
2881 : !! FUNCTION
2882 : !! Print info on the path.
2883 : !!
2884 : !! INPUTS
2885 : !! units=Unit numbers
2886 : !! [prtvol]=Verbosity level.
2887 : !! [header]=String to be printed as header for additional info.
2888 : !! [pre]=Optional string prepended to output e.g. #. Default: " "
2889 : !!
2890 : !! OUTPUT
2891 : !! Only printing
2892 : !!
2893 : !! SOURCE
2894 :
2895 20 : subroutine kpath_print(kpath, units, header,prtvol, pre)
2896 :
2897 : !Arguments ------------------------------------
2898 : !scalars
2899 : class(kpath_t),intent(in) :: kpath
2900 : integer,intent(in) :: units(:)
2901 : integer,optional,intent(in) :: prtvol
2902 : character(len=*),optional,intent(in) :: header,pre
2903 :
2904 : !Local variables-------------------------------
2905 : integer :: my_prtvol,ii
2906 : character(len=500) :: my_pre !, msg
2907 : ! *************************************************************************
2908 :
2909 10 : my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
2910 10 : my_pre = " "; if (present(pre)) my_pre = pre
2911 :
2912 10 : if (present(header)) call wrtout(units, sjoin(my_pre, '==== '//trim(adjustl(header))//' ==== '))
2913 :
2914 10 : call wrtout(units, sjoin(my_pre, " Number of points:", itoa(kpath%npts), ", ndivsmall:", itoa(kpath%ndivsm)))
2915 10 : call wrtout(units, sjoin(my_pre, " Boundaries and corresponding index in the k-points array:"))
2916 64 : do ii=1,kpath%nbounds
2917 64 : call wrtout(units, sjoin(my_pre, itoa(kpath%bounds2kpt(ii)), ktoa(kpath%bounds(:,ii))))
2918 : end do
2919 10 : call wrtout(units, " ")
2920 :
2921 10 : if (my_prtvol > 10) then
2922 0 : do ii=1,kpath%npts
2923 0 : call wrtout(units, sjoin(my_pre, ktoa(kpath%points(:,ii))))
2924 : end do
2925 : end if
2926 :
2927 10 : end subroutine kpath_print
2928 : !!***
2929 :
2930 : !----------------------------------------------------------------------
2931 :
2932 420599 : end module m_bz_mesh
2933 : !!***
|
