Branch data 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 : 1606 : subroutine kmesh_init(Kmesh, cryst, nkibz, kibz, kptopt, &
405 [ + + + - : 803 : 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 : 803 : integer,allocatable :: ktab(:),ktabi(:),ktabo(:)
424 : : real(dp) :: rm1t(3),kbz_wrap(3)
425 : 803 : real(dp),allocatable :: kbz(:,:),wtk(:)
426 : : ! *************************************************************************
427 : :
428 : : ! === Initial tests on input arguments ===
429 : 803 : ltest = (Cryst%timrev==1 .or. Cryst%timrev==2)
430 [ - + ]: 803 : ABI_CHECK(ltest, sjoin('Wrong value for timrev= ', itoa(Cryst%timrev)))
431 : :
432 [ + + ]: 803 : if (all(kptopt/= [1, 3])) then
433 : 14 : ABI_WARNING(sjoin("Not allowed value for kptopt: ", itoa(kptopt)))
434 : : end if
435 : :
436 : 803 : Kmesh%kptopt = kptopt
437 : 803 : nsym = Cryst%nsym; timrev = Cryst%timrev
438 : :
439 : : ! Find BZ from IBZ and fill tables ===
440 : 803 : nkbzX = nkibz*nsym*timrev ! Maximum possible number
441 [ - + + - : 2409 : ABI_MALLOC(kbz, (3,nkbzX))
+ - + - -
+ - + ]
442 [ - + + - : 2409 : ABI_MALLOC(wtk, (nkibz))
+ - + - -
+ - + ]
443 [ - + + - : 2409 : ABI_MALLOC(ktab, (nkbzX))
+ - + - -
+ - + ]
444 [ - + + - : 1606 : ABI_MALLOC(ktabi, (nkbzX))
- + - + ]
445 [ - + + - : 1606 : ABI_MALLOC(ktabo, (nkbzX))
- + - + ]
446 : :
447 [ + + ]: 803 : 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 : 781 : 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 [ + + ]: 803 : 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 : 803 : Kmesh%nbz = nkbz ! Number of points in the full BZ
471 : 803 : Kmesh%nibz = nkibz ! Number of points in the IBZ
472 : 803 : Kmesh%nsym = nsym ! Number of operations
473 : 803 : Kmesh%timrev= timrev ! 2 if time-reversal is used, 1 otherwise
474 : :
475 : : ! Arrays.
476 [ + + + + ]: 10439 : Kmesh%gmet = Cryst%gmet
477 [ + + + + ]: 10439 : Kmesh%gprimd = Cryst%gprimd
478 : :
479 [ - + + - : 2409 : ABI_MALLOC(Kmesh%bz ,(3,nkbz))
+ - + - -
+ - + -
+ ]
480 [ - + + + : 152790 : Kmesh%bz = kbz(:,1:nkbz ) ! Red. coordinates of points in full BZ.
+ + ]
481 [ - + + - : 2409 : ABI_MALLOC(Kmesh%ibz,(3,nkibz))
+ - + - -
+ - + -
+ ]
482 [ - + + + : 25010 : Kmesh%ibz = kibz(:,1:nkibz) ! Red. coordinates of points in IBZ.
+ + ]
483 : :
484 [ - + + - : 2409 : ABI_MALLOC(Kmesh%tab ,(nkbz))
+ - + - -
+ - + -
+ ]
485 [ - + + + ]: 39402 : Kmesh%tab = ktab (1:nkbz) ! Index of the irred. point in the array IBZ.
486 [ - + + - : 1606 : ABI_MALLOC(Kmesh%tabi,(nkbz))
- + - + -
+ ]
487 [ - + + + ]: 39402 : Kmesh%tabi = ktabi(1:nkbz) !-1 if time reversal must be used to obtain this point,
488 [ - + + - : 1606 : ABI_MALLOC(Kmesh%tabo,(nkbz))
- + - + -
+ ]
489 [ - + + + ]: 39402 : Kmesh%tabo = ktabo(1:nkbz) ! Symm. operation that rotates k_IBZ onto \pm k_BZ
490 : : ! (depending on tabi)
491 [ - + + - : 1606 : ABI_MALLOC(Kmesh%wt,(nkibz))
- + - + -
+ ]
492 [ + + ]: 6654 : Kmesh%wt(:)= wtk(1:nkibz) ! Weight for each k_IBZ
493 : :
494 [ - + + - : 4015 : ABI_MALLOC(Kmesh%rottbm1,(nkbz,timrev,nsym))
- + + - -
+ + - + -
+ - - + -
+ - + ]
495 [ - + + - : 3212 : ABI_MALLOC(Kmesh%rottb ,(nkbz,timrev,nsym))
- + + - -
+ + - - +
- + - + ]
496 : :
497 : 803 : do_hack = .FALSE.
498 [ + + ]: 803 : 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 : 781 : 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 [ - + + - : 2409 : ABI_MALLOC(Kmesh%umklp,(3,nkbz))
+ - + - -
+ - + -
+ ]
514 [ + + ]: 38599 : do ik_bz=1,nkbz
515 : 37796 : ik_ibz= Kmesh%tab (ik_bz)
516 : 37796 : isym = Kmesh%tabo(ik_bz)
517 : 37796 : itim = (3-Kmesh%tabi(ik_bz))/2
518 [ + + + + : 1248071 : 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 [ - + + - : 2409 : ABI_MALLOC(Kmesh%tabp,(nkbz))
+ - + - -
+ - + -
+ ]
522 [ + + ]: 38599 : do ik_bz=1,nkbz
523 : 37796 : isym =Kmesh%tabo(ik_bz)
524 : 37796 : ik_ibz=Kmesh%tab (ik_bz)
525 [ + + + + : 944900 : rm1t=MATMUL(TRANSPOSE(cryst%symrec(:,:,isym)),cryst%tnons(:,isym))
+ + + + +
- + - ]
526 [ + + ]: 151987 : Kmesh%tabp(ik_bz)=EXP(-(0.,1.)*two_pi*DOT_PRODUCT(kibz(:,ik_ibz),rm1t))
527 : : end do
528 : :
529 [ - + ]: 803 : ABI_FREE(kbz)
530 [ - + ]: 803 : ABI_FREE(wtk)
531 [ - + ]: 803 : ABI_FREE(ktab)
532 [ - + ]: 803 : ABI_FREE(ktabi)
533 [ - + ]: 803 : ABI_FREE(ktabo)
534 : :
535 : 803 : 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 : 781 : 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 : 1562 : integer :: g0(3), iperm(nbz),shlim(nbz+1)
715 : 1562 : 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 [ + + ]: 37449 : do ik=1,nbz
722 [ + + ]: 146672 : call wrap2_pmhalf(kbz(:,ik), kwrap, shift)
723 : 36668 : knorm(ik) = normv(kwrap,gmet,"G")
724 : 37449 : iperm(ik)= ik
725 : : end do
726 : :
727 : 781 : call sort_dp(nbz,knorm,iperm,KTOL)
728 : : !
729 : : ! The index of the initial (sorted) k-point in each shell
730 : 781 : nsh=1; norm_old=knorm(1)
731 : :
732 : 781 : shlim(1)=1; shlen(1) = norm_old
733 [ + + ]: 36668 : do ik_st=2,nbz
734 : 35887 : norm = knorm(ik_st)
735 [ + + ]: 36668 : if (ABS(norm-norm_old) > KTOL) then
736 : 5201 : norm_old = norm
737 : 5201 : nsh = nsh+1
738 : 5201 : shlim(nsh) = ik_st
739 : 5201 : shlen(nsh) = norm
740 : : end if
741 : : end do
742 : 781 : shlim(nsh+1)=nbz+1
743 : 781 : 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 : 781 : isok=.TRUE.
750 [ + + ]: 37449 : do ik=1,nbz
751 [ + + ]: 146672 : kbase(:)=kbz(:,ik)
752 : : !
753 [ + + ]: 110545 : do itim=1,timrev
754 [ + + ]: 3445188 : do isym=1,nsym
755 [ + + + + : 93391872 : krot(:)=(3-2*itim)*MATMUL(symrec(:,:,isym),kbase)
+ + + + +
+ ]
756 : :
757 : 206326624 : found=.FALSE.
758 : : #if 1
759 : : ! Old code
760 [ + - ]: 206326624 : do ikp=1,nbz
761 [ + + ]: 206326624 : if (isamek(krot,kbz(:,ikp),g0)) then
762 : 3335424 : found=.TRUE.
763 : 3335424 : krottb (ik ,itim,isym)=ikp
764 : 3335424 : 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 : 73096 : 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 [ - + ]: 781 : if (.not.isok) then
801 : 0 : ABI_ERROR('k-mesh not closed')
802 : : end if
803 : :
804 : : DBG_EXIT("COLL")
805 : :
806 : 781 : 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 : 465637433 : 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 [ + + ]: 1862549732 : isamek = isinteger(k1 - k2, TOL_KDIFF)
1022 : :
1023 [ + + ]: 465637433 : if (isamek) then
1024 [ + + ]: 19529808 : g0 = NINT(k1-k2)
1025 : : else
1026 [ + + ]: 1843019924 : g0 = HUGE(1)
1027 : : end if
1028 : :
1029 : 465637433 : 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 [ + + + - : 803 : 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 : 803 : is_irred_set=.TRUE.
1402 [ + + ]: 5851 : do ik1=1,nkibz-1
1403 [ + + ]: 20192 : k1=kibz(:,ik1)
1404 [ + + ]: 58632 : do ik2=ik1+1,nkibz
1405 [ + + ]: 211124 : k2=kibz(:,ik2)
1406 : :
1407 [ + + ]: 163299 : do itim=1,timrev
1408 [ + + ]: 5137571 : do isym=1,nsym
1409 [ + + ]: 4979320 : if (symafm(isym)==-1) CYCLE
1410 [ + + + + : 138598432 : knew = (3-2*itim) * MATMUL(symrec(:,:,isym),k2)
+ + + + +
+ ]
1411 [ - + ]: 5055414 : 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 [ - + ]: 803 : 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 : 803 : nkbz = 0
1431 [ + + ]: 6654 : do ikibz=1,nkibz
1432 : 5851 : 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 [ + + ]: 18288 : do itim=1,timrev
1437 [ + + ]: 538213 : do isym=1,nsym
1438 [ + + ]: 520728 : if (symafm(isym)==-1) CYCLE
1439 : : !
1440 : : ! * Form IS k
1441 [ + + + + : 14354592 : knew=(3-2*itim)*MATMUL(symrec(:,:,isym),kibz(:,ikibz))
+ + + + +
+ ]
1442 : : !
1443 : : ! * Check whether it has already been found (to within a RL vector).
1444 : 512664 : iold=0
1445 [ + + ]: 179734468 : do ikbz=1,nkbz
1446 [ + + ]: 179734468 : if (isamek(knew,kbz(:,ikbz),g0)) then
1447 : 512664 : 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 : 11634 : if (iold==0) then
1454 : 90812 : nkbz=nkbz+1
1455 : 90812 : wtk(ikibz)=wtk(ikibz)+one
1456 [ - + ]: 90812 : if (nkbz>nkbzmx) then
1457 : 0 : ABI_BUG(sjoin('nkbzmx too small, nkbzmx = ',itoa(nkbzmx),', increase nkbzmx !'))
1458 : : end if
1459 [ + + ]: 363248 : kbz(:,nkbz) = knew(:)
1460 : 90812 : ktab (nkbz) = ikibz
1461 : 90812 : ktabo(nkbz) = isym
1462 : 90812 : ktabi(nkbz) = 3-2*itim
1463 : : end if
1464 : : !
1465 : : end do
1466 : : end do
1467 : :
1468 : : end do !ikibz
1469 : :
1470 [ + + ]: 803 : 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 [ + + + + ]: 12505 : wtk = wtk/SUM(wtk)
1521 : :
1522 : : DBG_EXIT("COLL")
1523 : :
1524 : 803 : 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 [ - + - - : 420439 : end module m_bz_mesh
+ + - + +
- - + - +
+ - - + -
+ - + # #
# # # # #
# # # #
# ]
2933 : : !!***
|
