Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief K-points and crystal symmetry routines
10 : !> \author jgh
11 : ! **************************************************************************************************
12 : MODULE cryssym
13 :
14 : USE bibliography, ONLY: Togo2018,&
15 : Worlton1972,&
16 : cite_reference
17 : USE kinds, ONLY: dp
18 : USE kpsym, ONLY: group1s
19 : USE mathlib, ONLY: inv_3x3
20 : USE spglib_f08, ONLY: spg_get_international,&
21 : spg_get_major_version,&
22 : spg_get_micro_version,&
23 : spg_get_minor_version,&
24 : spg_get_multiplicity,&
25 : spg_get_pointgroup,&
26 : spg_get_schoenflies,&
27 : spg_get_symmetry
28 : USE string_utilities, ONLY: strip_control_codes
29 : #include "./base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 : PRIVATE
33 : PUBLIC :: csym_type, release_csym_type, print_crys_symmetry, print_kp_symmetry
34 : PUBLIC :: crys_sym_gen, kpoint_gen, kpoint_gen_general
35 :
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cryssym'
37 :
38 : ! **************************************************************************************************
39 : !> \brief CSM type
40 : !> \par Content:
41 : !>
42 : ! **************************************************************************************************
43 : TYPE csym_type
44 : LOGICAL :: symlib = .FALSE.
45 : LOGICAL :: fullgrid = .FALSE.
46 : LOGICAL :: inversion_only = .FALSE.
47 : LOGICAL :: spglib_reduction = .FALSE.
48 : LOGICAL :: spglib_backend = .FALSE.
49 : LOGICAL :: spglib_requested = .TRUE.
50 : INTEGER :: plevel = 0
51 : INTEGER :: punit = -1
52 : INTEGER :: istriz = -1
53 : REAL(KIND=dp) :: delta = 1.0e-8_dp
54 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp
55 : ! KPOINTS
56 : REAL(KIND=dp), DIMENSION(3) :: wvk0 = 0.0_dp
57 : INTEGER, DIMENSION(3) :: mesh = 0
58 : INTEGER :: nkpoint = 0
59 : INTEGER :: nat = 0
60 : INTEGER, DIMENSION(:), ALLOCATABLE :: atype
61 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: scoord
62 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: xkpoint
63 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: wkpoint
64 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: kpmesh
65 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: kplink
66 : INTEGER, DIMENSION(:), ALLOCATABLE :: kpop
67 : !SPGLIB
68 : CHARACTER(len=11) :: international_symbol = ""
69 : CHARACTER(len=6) :: pointgroup_symbol = ""
70 : CHARACTER(len=10) :: schoenflies = ""
71 : INTEGER :: n_operations = 0
72 : INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: rotations
73 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: translations
74 : !K290
75 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: rt
76 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: vt
77 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0
78 : INTEGER :: nrtot = 0
79 : INTEGER, DIMENSION(:), ALLOCATABLE :: ibrot
80 : END TYPE csym_type
81 :
82 : CONTAINS
83 :
84 : ! **************************************************************************************************
85 : !> \brief Release the CSYM type
86 : !> \param csym The CSYM type
87 : ! **************************************************************************************************
88 3746 : SUBROUTINE release_csym_type(csym)
89 : TYPE(csym_type) :: csym
90 :
91 3746 : IF (ALLOCATED(csym%rotations)) THEN
92 3463 : DEALLOCATE (csym%rotations)
93 : END IF
94 3746 : IF (ALLOCATED(csym%translations)) THEN
95 3463 : DEALLOCATE (csym%translations)
96 : END IF
97 3746 : IF (ALLOCATED(csym%atype)) THEN
98 3746 : DEALLOCATE (csym%atype)
99 : END IF
100 3746 : IF (ALLOCATED(csym%scoord)) THEN
101 3746 : DEALLOCATE (csym%scoord)
102 : END IF
103 3746 : IF (ALLOCATED(csym%xkpoint)) THEN
104 3456 : DEALLOCATE (csym%xkpoint)
105 : END IF
106 3746 : IF (ALLOCATED(csym%wkpoint)) THEN
107 3456 : DEALLOCATE (csym%wkpoint)
108 : END IF
109 3746 : IF (ALLOCATED(csym%kpmesh)) THEN
110 3456 : DEALLOCATE (csym%kpmesh)
111 : END IF
112 3746 : IF (ALLOCATED(csym%kplink)) THEN
113 3456 : DEALLOCATE (csym%kplink)
114 : END IF
115 3746 : IF (ALLOCATED(csym%kpop)) THEN
116 3456 : DEALLOCATE (csym%kpop)
117 : END IF
118 3746 : IF (ALLOCATED(csym%rt)) THEN
119 3456 : DEALLOCATE (csym%rt)
120 : END IF
121 3746 : IF (ALLOCATED(csym%vt)) THEN
122 3456 : DEALLOCATE (csym%vt)
123 : END IF
124 3746 : IF (ALLOCATED(csym%f0)) THEN
125 3456 : DEALLOCATE (csym%f0)
126 : END IF
127 3746 : IF (ALLOCATED(csym%ibrot)) THEN
128 3456 : DEALLOCATE (csym%ibrot)
129 : END IF
130 :
131 3746 : END SUBROUTINE release_csym_type
132 :
133 : ! **************************************************************************************************
134 : !> \brief ...
135 : !> \param csym ...
136 : !> \param scoor ...
137 : !> \param types ...
138 : !> \param hmat ...
139 : !> \param delta ...
140 : !> \param iounit ...
141 : !> \param use_spglib ...
142 : ! **************************************************************************************************
143 3746 : SUBROUTINE crys_sym_gen(csym, scoor, types, hmat, delta, iounit, use_spglib)
144 : TYPE(csym_type) :: csym
145 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scoor
146 : INTEGER, DIMENSION(:), INTENT(IN) :: types
147 : REAL(KIND=dp), INTENT(IN) :: hmat(3, 3)
148 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: delta
149 : INTEGER, INTENT(IN), OPTIONAL :: iounit
150 : LOGICAL, INTENT(IN), OPTIONAL :: use_spglib
151 :
152 : CHARACTER(LEN=*), PARAMETER :: routineN = 'crys_sym_gen'
153 :
154 : INTEGER :: handle, ierr, major, micro, minor, nat, &
155 : nop, tra_mat(3, 3)
156 : LOGICAL :: my_use_spglib, spglib
157 :
158 3746 : CALL timeset(routineN, handle)
159 :
160 : !..total number of atoms
161 3746 : nat = SIZE(scoor, 2)
162 3746 : csym%nat = nat
163 :
164 : ! output unit
165 3746 : IF (PRESENT(iounit)) THEN
166 3734 : csym%punit = iounit
167 : ELSE
168 12 : csym%punit = -1
169 : END IF
170 :
171 : ! accuracy for symmetry
172 3746 : IF (PRESENT(delta)) THEN
173 3746 : csym%delta = delta
174 : ELSE
175 0 : csym%delta = 1.e-6_dp
176 : END IF
177 :
178 : !..set cell values
179 48698 : csym%hmat = hmat
180 :
181 : ! atom types
182 11238 : ALLOCATE (csym%atype(nat))
183 24985 : csym%atype(1:nat) = types(1:nat)
184 :
185 : ! scaled coordinates
186 11238 : ALLOCATE (csym%scoord(3, nat))
187 88702 : csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
188 :
189 3746 : csym%n_operations = 0
190 :
191 : !..try spglib
192 3746 : my_use_spglib = .TRUE.
193 3746 : IF (PRESENT(use_spglib)) my_use_spglib = use_spglib
194 290 : csym%spglib_requested = my_use_spglib
195 3456 : IF (.NOT. my_use_spglib) THEN
196 : spglib = .FALSE.
197 : ELSE
198 3464 : major = spg_get_major_version()
199 3464 : minor = spg_get_minor_version()
200 3464 : micro = spg_get_micro_version()
201 3464 : IF (major == 0) THEN
202 0 : CALL cp_warn(__LOCATION__, "Symmetry library SPGLIB not available")
203 0 : spglib = .FALSE.
204 : ELSE
205 3464 : spglib = .TRUE.
206 3464 : CALL cite_reference(Togo2018)
207 3464 : ierr = spg_get_international(csym%international_symbol, TRANSPOSE(hmat), scoor, types, nat, delta)
208 3464 : IF (ierr == 0) THEN
209 1 : CALL cp_warn(__LOCATION__, "Symmetry Library SPGLIB failed")
210 1 : spglib = .FALSE.
211 : ELSE
212 3463 : nop = spg_get_multiplicity(TRANSPOSE(hmat), scoor, types, nat, delta)
213 17315 : ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
214 3463 : csym%n_operations = nop
215 : ierr = spg_get_symmetry(csym%rotations, csym%translations, nop, &
216 3463 : TRANSPOSE(hmat), scoor, types, nat, delta)
217 : ! Schoenflies Symbol
218 3463 : csym%schoenflies = ' '
219 3463 : ierr = spg_get_schoenflies(csym%schoenflies(1:7), TRANSPOSE(hmat), scoor, types, nat, delta)
220 : ! Point Group
221 3463 : csym%pointgroup_symbol = ' '
222 3463 : tra_mat = 0
223 : ierr = spg_get_pointgroup(csym%pointgroup_symbol, tra_mat, &
224 3463 : csym%rotations, csym%n_operations)
225 :
226 3463 : CALL strip_control_codes(csym%international_symbol)
227 3463 : CALL strip_control_codes(csym%schoenflies)
228 3463 : CALL strip_control_codes(csym%pointgroup_symbol)
229 : END IF
230 : END IF
231 : END IF
232 3746 : csym%symlib = spglib
233 :
234 3746 : CALL timestop(handle)
235 :
236 3746 : END SUBROUTINE crys_sym_gen
237 :
238 : ! **************************************************************************************************
239 : !> \brief ...
240 : !> \param csym ...
241 : !> \param nk ...
242 : !> \param symm ...
243 : !> \param shift ...
244 : !> \param full_grid ...
245 : !> \param gamma_centered ...
246 : !> \param inversion_symmetry_only ...
247 : !> \param use_spglib_reduction ...
248 : !> \param use_spglib_backend ...
249 : ! **************************************************************************************************
250 3430 : SUBROUTINE kpoint_gen(csym, nk, symm, shift, full_grid, gamma_centered, &
251 : inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
252 : TYPE(csym_type) :: csym
253 : INTEGER, INTENT(IN) :: nk(3)
254 : LOGICAL, INTENT(IN), OPTIONAL :: symm
255 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: shift(3)
256 : LOGICAL, INTENT(IN), OPTIONAL :: full_grid, gamma_centered, &
257 : inversion_symmetry_only, &
258 : use_spglib_reduction, &
259 : use_spglib_backend
260 :
261 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_gen'
262 :
263 : INTEGER :: handle, i, ik, j, nkp, nkpts
264 3430 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kpop, rep_to_red
265 : LOGICAL :: fullmesh, gamma_mesh, inversion_only, &
266 : spglib_backend, spglib_reduction
267 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkp
268 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
269 :
270 3430 : CALL timeset(routineN, handle)
271 :
272 3430 : IF (PRESENT(shift)) THEN
273 13720 : csym%wvk0 = shift
274 : ELSE
275 0 : csym%wvk0 = 0.0_dp
276 : END IF
277 :
278 3430 : csym%istriz = -1
279 3430 : IF (PRESENT(symm)) THEN
280 3430 : IF (symm) csym%istriz = 1
281 : END IF
282 :
283 3430 : IF (PRESENT(full_grid)) THEN
284 3418 : fullmesh = full_grid
285 : ELSE
286 : fullmesh = .FALSE.
287 : END IF
288 3430 : csym%fullgrid = fullmesh
289 :
290 3430 : IF (PRESENT(gamma_centered)) THEN
291 3418 : gamma_mesh = gamma_centered
292 : ELSE
293 12 : gamma_mesh = .FALSE.
294 : END IF
295 :
296 3430 : IF (PRESENT(inversion_symmetry_only)) THEN
297 3418 : inversion_only = inversion_symmetry_only
298 : ELSE
299 : inversion_only = .FALSE.
300 : END IF
301 3430 : csym%inversion_only = inversion_only
302 :
303 3430 : IF (PRESENT(use_spglib_reduction)) THEN
304 3418 : spglib_reduction = use_spglib_reduction
305 : ELSE
306 12 : spglib_reduction = .FALSE.
307 : END IF
308 3430 : csym%spglib_reduction = spglib_reduction
309 :
310 3430 : IF (PRESENT(use_spglib_backend)) THEN
311 3418 : spglib_backend = use_spglib_backend
312 : ELSE
313 : spglib_backend = .FALSE.
314 : END IF
315 3430 : csym%spglib_backend = spglib_backend
316 :
317 3430 : IF (spglib_backend .AND. .NOT. spglib_reduction) THEN
318 : CALL cp_abort(__LOCATION__, &
319 0 : "SYMMETRY_BACKEND SPGLIB requires SYMMETRY_REDUCTION_METHOD SPGLIB")
320 : END IF
321 : IF (csym%istriz == 1 .AND. .NOT. fullmesh .AND. .NOT. inversion_only .AND. &
322 3430 : (spglib_backend .OR. spglib_reduction) .AND. .NOT. csym%symlib) THEN
323 : CALL cp_abort(__LOCATION__, &
324 0 : "SPGLIB k-point symmetry was requested, but SPGLIB is not available")
325 : END IF
326 :
327 3430 : csym%nkpoint = 0
328 13720 : csym%mesh(1:3) = nk(1:3)
329 3430 : csym%nrtot = 0
330 3430 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
331 3430 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
332 3430 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
333 3430 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
334 6860 : ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
335 :
336 3430 : nkpts = nk(1)*nk(2)*nk(3)
337 24010 : ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
338 : ! kp: link
339 10290 : ALLOCATE (csym%kplink(2, nkpts))
340 106834 : csym%kplink = 0
341 3430 : kpop = 0
342 :
343 : ! go through all the options
344 3430 : IF (csym%symlib) THEN
345 : ! symmetry library is available
346 3152 : IF (fullmesh) THEN
347 : ! full mesh requested
348 346 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
349 346 : IF (csym%istriz == 1) THEN
350 : ! use inversion symmetry
351 346 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
352 : ELSE
353 : ! full kpoint mesh is used
354 : END IF
355 2806 : ELSE IF (csym%istriz /= 1 .OR. inversion_only) THEN
356 : ! use inversion symmetry
357 1108 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
358 1108 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
359 : ELSE
360 : ! use symmetry library to reduce k-points
361 1698 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
362 1698 : IF (spglib_backend) THEN
363 760 : CALL kp_symmetry_spglib(csym, xkp, wkp, kpop)
364 : ELSE
365 938 : CALL kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction=spglib_reduction)
366 : END IF
367 :
368 : END IF
369 : ELSE
370 : ! no symmetry library is available
371 278 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
372 278 : IF (csym%istriz == 1 .AND. .NOT. fullmesh .AND. .NOT. inversion_only) THEN
373 : ! fall back to the K290 atom mapping when SPGLIB is not linked
374 12 : CALL kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction=.FALSE.)
375 266 : ELSE IF (csym%istriz /= 1 .AND. fullmesh) THEN
376 : ! full kpoint mesh is used
377 1148 : DO i = 1, nkpts
378 1148 : csym%kplink(1, i) = i
379 : END DO
380 : ELSE
381 : ! use inversion symmetry
382 164 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
383 : END IF
384 : END IF
385 : ! count kpoints
386 3430 : nkp = 0
387 37898 : DO i = 1, nkpts
388 37898 : IF (wkp(i) > 0.0_dp) nkp = nkp + 1
389 : END DO
390 :
391 : ! store reduced kpoint set
392 3430 : csym%nkpoint = nkp
393 17150 : ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
394 6860 : ALLOCATE (rep_to_red(nkpts))
395 3430 : rep_to_red(:) = 0
396 3430 : j = 0
397 37898 : DO ik = 1, nkpts
398 37898 : IF (wkp(ik) > 0.0_dp) THEN
399 12316 : j = j + 1
400 12316 : csym%wkpoint(j) = wkp(ik)
401 49264 : csym%xkpoint(1:3, j) = xkp(1:3, ik)
402 12316 : rep_to_red(ik) = j
403 : END IF
404 : END DO
405 3430 : CPASSERT(j == nkp)
406 :
407 : ! kp: mesh
408 6860 : ALLOCATE (csym%kpmesh(3, nkpts))
409 141302 : csym%kpmesh(1:3, 1:nkpts) = xkp(1:3, 1:nkpts)
410 :
411 : ! kp: link
412 37898 : DO ik = 1, nkpts
413 34468 : i = csym%kplink(1, ik)
414 34468 : CPASSERT(i > 0 .AND. i <= nkpts)
415 34468 : csym%kplink(2, ik) = rep_to_red(i)
416 37898 : CPASSERT(csym%kplink(2, ik) > 0)
417 : END DO
418 3430 : DEALLOCATE (rep_to_red)
419 :
420 : ! kp: operations
421 6860 : ALLOCATE (csym%kpop(nkpts))
422 3430 : IF (csym%nrtot > 0 .AND. csym%istriz == 1 .AND. .NOT. fullmesh .AND. &
423 : .NOT. inversion_only) THEN
424 : ! atomic symmetry operations possible
425 21506 : csym%kpop(1:nkpts) = kpop(1:nkpts)
426 21506 : DO ik = 1, nkpts
427 21506 : CPASSERT(csym%kpop(ik) /= 0)
428 : END DO
429 : ELSE
430 : ! only time reversal symmetry
431 16392 : DO ik = 1, nkpts
432 16392 : IF (wkp(ik) > 0.0_dp) THEN
433 8146 : csym%kpop(ik) = 1
434 : ELSE
435 6524 : csym%kpop(ik) = 2
436 : END IF
437 : END DO
438 : END IF
439 :
440 3430 : DEALLOCATE (xkp, wkp, kpop)
441 :
442 3430 : CALL timestop(handle)
443 :
444 3430 : END SUBROUTINE kpoint_gen
445 :
446 : ! **************************************************************************************************
447 : !> \brief Reduce an explicitly supplied GENERAL k-point set.
448 : !> \param csym ...
449 : !> \param xkp_in explicit k-point coordinates in reciprocal lattice coordinates
450 : !> \param wkp_in explicit k-point weights
451 : !> \param symm ...
452 : !> \param full_grid ...
453 : !> \param inversion_symmetry_only ...
454 : !> \param use_spglib_reduction ...
455 : !> \param use_spglib_backend ...
456 : ! **************************************************************************************************
457 26 : SUBROUTINE kpoint_gen_general(csym, xkp_in, wkp_in, symm, full_grid, &
458 : inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
459 : TYPE(csym_type) :: csym
460 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_in
461 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_in
462 : LOGICAL, INTENT(IN), OPTIONAL :: symm, full_grid, &
463 : inversion_symmetry_only, &
464 : use_spglib_reduction, &
465 : use_spglib_backend
466 :
467 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_gen_general'
468 :
469 : INTEGER :: handle, i, nfull
470 : LOGICAL :: atomic_symmetry, fullmesh, &
471 : inversion_only, spglib_backend, &
472 : spglib_reduction
473 : REAL(KIND=dp) :: weight_eps
474 :
475 26 : CALL timeset(routineN, handle)
476 :
477 26 : nfull = SIZE(wkp_in)
478 26 : CPASSERT(SIZE(xkp_in, 1) == 3)
479 26 : CPASSERT(SIZE(xkp_in, 2) == nfull)
480 :
481 26 : atomic_symmetry = .FALSE.
482 26 : IF (PRESENT(symm)) atomic_symmetry = symm
483 0 : csym%istriz = -1
484 26 : IF (atomic_symmetry) csym%istriz = 1
485 26 : fullmesh = .FALSE.
486 26 : IF (PRESENT(full_grid)) fullmesh = full_grid
487 26 : inversion_only = .FALSE.
488 26 : IF (PRESENT(inversion_symmetry_only)) inversion_only = inversion_symmetry_only
489 26 : spglib_reduction = .FALSE.
490 26 : IF (PRESENT(use_spglib_reduction)) spglib_reduction = use_spglib_reduction
491 26 : spglib_backend = .FALSE.
492 26 : IF (PRESENT(use_spglib_backend)) spglib_backend = use_spglib_backend
493 :
494 26 : csym%fullgrid = fullmesh
495 26 : csym%inversion_only = inversion_only
496 26 : csym%spglib_reduction = spglib_reduction
497 26 : csym%spglib_backend = spglib_backend
498 26 : csym%nkpoint = 0
499 104 : csym%mesh(1:3) = 0
500 26 : csym%nrtot = 0
501 26 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
502 26 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
503 26 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
504 26 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
505 52 : ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
506 26 : IF (ALLOCATED(csym%xkpoint)) DEALLOCATE (csym%xkpoint)
507 26 : IF (ALLOCATED(csym%wkpoint)) DEALLOCATE (csym%wkpoint)
508 26 : IF (ALLOCATED(csym%kpmesh)) DEALLOCATE (csym%kpmesh)
509 26 : IF (ALLOCATED(csym%kplink)) DEALLOCATE (csym%kplink)
510 26 : IF (ALLOCATED(csym%kpop)) DEALLOCATE (csym%kpop)
511 :
512 182 : ALLOCATE (csym%kpmesh(3, nfull), csym%kplink(2, nfull), csym%kpop(nfull))
513 858 : csym%kpmesh(1:3, 1:nfull) = xkp_in(1:3, 1:nfull)
514 650 : csym%kplink = 0
515 234 : csym%kpop = 1
516 :
517 26 : IF (.NOT. atomic_symmetry .OR. fullmesh) THEN
518 0 : csym%nkpoint = nfull
519 0 : ALLOCATE (csym%xkpoint(3, nfull), csym%wkpoint(nfull))
520 0 : csym%xkpoint(1:3, 1:nfull) = xkp_in(1:3, 1:nfull)
521 0 : csym%wkpoint(1:nfull) = wkp_in(1:nfull)
522 0 : DO i = 1, nfull
523 0 : csym%kplink(1:2, i) = i
524 : END DO
525 26 : ELSE IF (inversion_only) THEN
526 0 : CALL reduce_general_inversion(csym, xkp_in, wkp_in)
527 : ELSE
528 26 : weight_eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
529 234 : IF (ANY(ABS(wkp_in(1:nfull) - wkp_in(1)) > weight_eps)) THEN
530 : CALL cp_abort(__LOCATION__, &
531 0 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires equal explicit weights.")
532 : END IF
533 26 : IF (spglib_backend) THEN
534 18 : IF (.NOT. csym%symlib) THEN
535 : CALL cp_abort(__LOCATION__, &
536 0 : "SCHEME GENERAL with SYMMETRY_BACKEND SPGLIB requires SPGLIB.")
537 : END IF
538 18 : CALL reduce_general_spglib(csym, xkp_in)
539 8 : ELSE IF (spglib_reduction) THEN
540 4 : IF (.NOT. csym%symlib) THEN
541 : CALL cp_abort(__LOCATION__, &
542 0 : "SCHEME GENERAL with SYMMETRY_REDUCTION_METHOD SPGLIB requires SPGLIB.")
543 : END IF
544 4 : CALL setup_k290_operations(csym)
545 4 : CALL reduce_general_spglib_k290(csym, xkp_in)
546 : ELSE
547 4 : CALL setup_k290_operations(csym)
548 4 : CALL reduce_general_k290(csym, xkp_in)
549 : END IF
550 : END IF
551 :
552 26 : CALL timestop(handle)
553 :
554 26 : END SUBROUTINE kpoint_gen_general
555 :
556 : ! **************************************************************************************************
557 : !> \brief ...
558 : !> \param csym ...
559 : !> \param xkp ...
560 : !> \param wkp ...
561 : !> \param kpop ...
562 : !> \param use_spglib_reduction ...
563 : ! **************************************************************************************************
564 950 : SUBROUTINE kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction)
565 : TYPE(csym_type) :: csym
566 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
567 : REAL(KIND=dp), DIMENSION(:) :: wkp
568 : INTEGER, DIMENSION(:) :: kpop
569 : LOGICAL, INTENT(IN), OPTIONAL :: use_spglib_reduction
570 :
571 : INTEGER :: nc, nrot
572 950 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
573 : INTEGER, DIMENSION(48) :: ib
574 : LOGICAL :: spglib_reduction
575 : REAL(KIND=dp) :: alat
576 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3
577 : REAL(KIND=dp), DIMENSION(3, 3, 48) :: r
578 :
579 950 : IF (PRESENT(use_spglib_reduction)) THEN
580 950 : spglib_reduction = use_spglib_reduction
581 : ELSE
582 : spglib_reduction = .FALSE.
583 : END IF
584 :
585 950 : CALL detect_k290_symmetry(csym, a1, a2, a3, b1, b2, b3, alat, nc, ib, r)
586 950 : IF (spglib_reduction) CALL cite_reference(Togo2018)
587 :
588 950 : IF (csym%n_operations > nc .AND. .NOT. spglib_reduction) THEN
589 702 : ALLOCATE (srot(3, 3, csym%n_operations))
590 234 : CALL setup_spglib_operations(csym, srot, nrot)
591 234 : CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
592 468 : DEALLOCATE (srot)
593 708 : ELSE IF (spglib_reduction) THEN
594 48 : ALLOCATE (srot(3, 3, csym%n_operations))
595 16 : CALL setup_spglib_reduction_rotations(csym, srot, nrot)
596 : CALL reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
597 16 : a1, a2, a3, b1, b2, b3, alat)
598 32 : DEALLOCATE (srot)
599 : ELSE
600 700 : CALL reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
601 : END IF
602 :
603 950 : END SUBROUTINE kp_symmetry
604 :
605 : ! **************************************************************************************************
606 : !> \brief Reduce a CP2K Monkhorst-Pack mesh using SPGLIB symmetry operations
607 : !> \param csym ...
608 : !> \param xkp ...
609 : !> \param wkp ...
610 : !> \param kpop ...
611 : ! **************************************************************************************************
612 760 : SUBROUTINE kp_symmetry_spglib(csym, xkp, wkp, kpop)
613 : TYPE(csym_type) :: csym
614 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
615 : REAL(KIND=dp), DIMENSION(:) :: wkp
616 : INTEGER, DIMENSION(:) :: kpop
617 :
618 : INTEGER :: iou, nrot
619 760 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
620 :
621 760 : iou = csym%punit
622 760 : IF (iou > 0) THEN
623 : WRITE (iou, '(/,(T2,A79))') &
624 321 : "*******************************************************************************", &
625 321 : "** Special K-Point Generation by SPGLIB **", &
626 642 : "*******************************************************************************"
627 : END IF
628 760 : CALL cite_reference(Togo2018)
629 :
630 2280 : ALLOCATE (srot(3, 3, csym%n_operations))
631 760 : CALL setup_spglib_operations(csym, srot, nrot)
632 760 : CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
633 760 : DEALLOCATE (srot)
634 :
635 760 : IF (iou > 0) THEN
636 : WRITE (iou, '((T2,A79))') &
637 321 : "*******************************************************************************", &
638 321 : "** Finished SPGLIB **", &
639 642 : "*******************************************************************************"
640 : END IF
641 :
642 760 : END SUBROUTINE kp_symmetry_spglib
643 :
644 : ! **************************************************************************************************
645 : !> \brief Store K290 atomic symmetry operations.
646 : !> \param csym ...
647 : ! **************************************************************************************************
648 8 : SUBROUTINE setup_k290_operations(csym)
649 : TYPE(csym_type) :: csym
650 :
651 : INTEGER :: nc
652 : INTEGER, DIMENSION(48) :: ib
653 : REAL(KIND=dp) :: alat
654 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3
655 : REAL(KIND=dp), DIMENSION(3, 3, 48) :: r
656 :
657 8 : CALL detect_k290_symmetry(csym, a1, a2, a3, b1, b2, b3, alat, nc, ib, r)
658 :
659 8 : END SUBROUTINE setup_k290_operations
660 :
661 : ! **************************************************************************************************
662 : !> \brief Detect and store atomic symmetry operations with the K290/Worlton routines.
663 : !> \param csym crystal symmetry data
664 : !> \param a1 first lattice vector
665 : !> \param a2 second lattice vector
666 : !> \param a3 third lattice vector
667 : !> \param b1 first reciprocal lattice vector
668 : !> \param b2 second reciprocal lattice vector
669 : !> \param b3 third reciprocal lattice vector
670 : !> \param alat lattice scaling used by K290
671 : !> \param nc number of detected operations
672 : !> \param ib indices of the detected rotations
673 : !> \param r K290 rotation matrices
674 : ! **************************************************************************************************
675 958 : SUBROUTINE detect_k290_symmetry(csym, a1, a2, a3, b1, b2, b3, alat, nc, ib, r)
676 : TYPE(csym_type), INTENT(INOUT) :: csym
677 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: a1, a2, a3, b1, b2, b3
678 : REAL(KIND=dp), INTENT(OUT) :: alat
679 : INTEGER, INTENT(OUT) :: nc
680 : INTEGER, DIMENSION(48), INTENT(OUT) :: ib
681 : REAL(KIND=dp), DIMENSION(3, 3, 48), INTENT(OUT) :: r
682 :
683 : INTEGER :: i, ihc, ihg, indpg, iou, isy, li, nat, &
684 : ntvec
685 958 : INTEGER, ALLOCATABLE, DIMENSION(:) :: isc, ty
686 958 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0
687 958 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rx, tvec, xkapa
688 : REAL(KIND=dp), DIMENSION(3) :: origin
689 : REAL(KIND=dp), DIMENSION(3, 48) :: vt
690 :
691 958 : iou = csym%punit
692 958 : nat = csym%nat
693 958 : CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
694 :
695 9580 : ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, nat), ty(nat), isc(nat), f0(49, nat))
696 6050 : ty(1:nat) = csym%atype(1:nat)
697 6050 : DO i = 1, nat
698 87522 : xkapa(1:3, i) = MATMUL(csym%hmat, csym%scoord(1:3, i))
699 : END DO
700 :
701 958 : IF (iou > 0) THEN
702 : WRITE (iou, '(/,(T2,A79))') &
703 335 : "*******************************************************************************", &
704 335 : "** Crystal Symmetry Analysis by K290 **", &
705 670 : "*******************************************************************************"
706 : END IF
707 958 : CALL cite_reference(Worlton1972)
708 :
709 : CALL GROUP1s(iou, a1, a2, a3, nat, ty, xkapa, b1, b2, b3, &
710 : ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
711 958 : vt, f0, r, tvec, origin, rx, isc, csym%delta)
712 :
713 958 : IF (iou > 0) THEN
714 : WRITE (iou, '((T2,A79))') &
715 335 : "*******************************************************************************", &
716 335 : "** Finished K290 **", &
717 670 : "*******************************************************************************"
718 : END IF
719 :
720 958 : csym%nrtot = nc
721 958 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
722 958 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
723 958 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
724 958 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
725 9580 : ALLOCATE (csym%rt(3, 3, nc), csym%vt(3, nc), csym%ibrot(nc), csym%f0(nat, nc))
726 65934 : csym%vt(1:3, 1:nc) = vt(1:3, 1:nc)
727 17202 : DO i = 1, nc
728 211172 : csym%rt(1:3, 1:3, i) = r(1:3, 1:3, ib(i))
729 92034 : csym%f0(1:nat, i) = f0(i, 1:nat)
730 : END DO
731 17202 : csym%ibrot(1:nc) = ib(1:nc)
732 :
733 958 : DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
734 :
735 958 : END SUBROUTINE detect_k290_symmetry
736 :
737 : ! **************************************************************************************************
738 : !> \brief Return K290 lattice vectors and reciprocal vectors.
739 : !> \param csym ...
740 : !> \param a1 first lattice vector
741 : !> \param a2 second lattice vector
742 : !> \param a3 third lattice vector
743 : !> \param b1 first reciprocal lattice vector
744 : !> \param b2 second reciprocal lattice vector
745 : !> \param b3 third reciprocal lattice vector
746 : !> \param alat lattice scaling used by K290
747 : ! **************************************************************************************************
748 966 : SUBROUTINE setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
749 : TYPE(csym_type), INTENT(IN) :: csym
750 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: a1, a2, a3, b1, b2, b3
751 : REAL(KIND=dp), INTENT(OUT) :: alat
752 :
753 : REAL(KIND=dp) :: volum
754 :
755 3864 : a1(1:3) = csym%hmat(1:3, 1)
756 3864 : a2(1:3) = csym%hmat(1:3, 2)
757 3864 : a3(1:3) = csym%hmat(1:3, 3)
758 3864 : alat = SQRT(SUM(a1**2))
759 : volum = a1(1)*a2(2)*a3(3) + a2(1)*a3(2)*a1(3) + &
760 : a3(1)*a1(2)*a2(3) - a1(3)*a2(2)*a3(1) - &
761 966 : a2(3)*a3(2)*a1(1) - a3(3)*a1(2)*a2(1)
762 966 : volum = ABS(volum)
763 966 : b1(1) = (a2(2)*a3(3) - a2(3)*a3(2))/volum
764 966 : b1(2) = (a2(3)*a3(1) - a2(1)*a3(3))/volum
765 966 : b1(3) = (a2(1)*a3(2) - a2(2)*a3(1))/volum
766 966 : b2(1) = (a3(2)*a1(3) - a3(3)*a1(2))/volum
767 966 : b2(2) = (a3(3)*a1(1) - a3(1)*a1(3))/volum
768 966 : b2(3) = (a3(1)*a1(2) - a3(2)*a1(1))/volum
769 966 : b3(1) = (a1(2)*a2(3) - a1(3)*a2(2))/volum
770 966 : b3(2) = (a1(3)*a2(1) - a1(1)*a2(3))/volum
771 966 : b3(3) = (a1(1)*a2(2) - a1(2)*a2(1))/volum
772 :
773 966 : END SUBROUTINE setup_k290_lattice
774 :
775 : ! **************************************************************************************************
776 : !> \brief Store usable SPGLIB space-group operations for k-point symmetry
777 : !> \param csym ...
778 : !> \param srot integer rotations in fractional coordinates
779 : !> \param nrot number of stored rotations
780 : ! **************************************************************************************************
781 1012 : SUBROUTINE setup_spglib_operations(csym, srot, nrot)
782 : TYPE(csym_type) :: csym
783 : INTEGER, DIMENSION(:, :, :), INTENT(OUT) :: srot
784 : INTEGER, INTENT(OUT) :: nrot
785 :
786 : INTEGER :: iop, jop, pass
787 1012 : INTEGER, ALLOCATABLE, DIMENSION(:) :: perm
788 : INTEGER, DIMENSION(3, 3) :: eye, frot, irot
789 : LOGICAL :: duplicate, identity, valid, &
790 : zero_translation
791 : REAL(KIND=dp) :: eps
792 : REAL(KIND=dp), DIMENSION(3, 3) :: h_inv, rfrac
793 :
794 1012 : CPASSERT(csym%symlib)
795 :
796 652156 : srot = 0
797 1012 : csym%nrtot = 0
798 1012 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
799 1012 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
800 1012 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
801 1012 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
802 5060 : ALLOCATE (csym%rt(3, 3, csym%n_operations), csym%vt(3, csym%n_operations))
803 6072 : ALLOCATE (csym%ibrot(csym%n_operations), csym%f0(csym%nat, csym%n_operations))
804 652156 : csym%rt = 0.0_dp
805 201364 : csym%vt = 0.0_dp
806 51100 : csym%ibrot = 0
807 417604 : csym%f0 = 0
808 :
809 1012 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
810 1012 : h_inv = inv_3x3(csym%hmat)
811 3036 : ALLOCATE (perm(csym%nat))
812 :
813 1012 : eye = 0
814 1012 : eye(1, 1) = 1
815 1012 : eye(2, 2) = 1
816 1012 : eye(3, 3) = 1
817 :
818 1012 : nrot = 0
819 : ! Operation 1 is used as the untransformed representative k-point.
820 : ! Prefer integer translations before fractional alternatives with the same rotation.
821 5060 : DO pass = 1, 4
822 205412 : DO iop = 1, csym%n_operations
823 2604576 : irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
824 2604576 : frot(1:3, 1:3) = TRANSPOSE(irot(1:3, 1:3))
825 458128 : identity = ALL(frot == eye)
826 : zero_translation = ALL(ABS(csym%translations(1:3, iop) - &
827 346456 : ANINT(csym%translations(1:3, iop))) <= eps)
828 200352 : IF (pass == 1 .AND. (.NOT. identity .OR. .NOT. zero_translation)) CYCLE
829 151276 : IF (pass == 2 .AND. (identity .OR. .NOT. zero_translation)) CYCLE
830 110490 : IF (pass == 3 .AND. (.NOT. identity .OR. zero_translation)) CYCLE
831 61202 : IF (pass == 4 .AND. (identity .OR. zero_translation)) CYCLE
832 :
833 50088 : duplicate = .FALSE.
834 1141146 : DO jop = 1, nrot
835 2786168 : IF (ALL(frot == srot(:, :, jop))) THEN
836 : duplicate = .TRUE.
837 : EXIT
838 : END IF
839 : END DO
840 50088 : IF (duplicate) CYCLE
841 :
842 16776 : CALL spglib_atom_permutation(csym, frot, csym%translations(:, iop), perm, valid)
843 16776 : IF (.NOT. valid) CYCLE
844 :
845 16776 : nrot = nrot + 1
846 :
847 218088 : srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
848 218088 : rfrac(1:3, 1:3) = REAL(frot(1:3, 1:3), KIND=dp)
849 1526616 : csym%rt(1:3, 1:3, nrot) = MATMUL(csym%hmat, MATMUL(rfrac, h_inv))
850 67104 : csym%vt(1:3, nrot) = csym%translations(1:3, iop)
851 16776 : csym%ibrot(nrot) = nrot
852 123136 : csym%f0(1:csym%nat, nrot) = perm(1:csym%nat)
853 : END DO
854 : END DO
855 :
856 1012 : DEALLOCATE (perm)
857 1012 : csym%nrtot = nrot
858 1012 : IF (nrot == 0) CALL cp_abort(__LOCATION__, "SPGLIB did not return usable symmetry operations")
859 :
860 1012 : END SUBROUTINE setup_spglib_operations
861 :
862 : ! **************************************************************************************************
863 : !> \brief Store unique SPGLIB rotations for K290-backend diagnostic reduction
864 : !> \param csym ...
865 : !> \param srot integer rotations in fractional coordinates
866 : !> \param nrot number of stored rotations
867 : ! **************************************************************************************************
868 20 : SUBROUTINE setup_spglib_reduction_rotations(csym, srot, nrot)
869 : TYPE(csym_type) :: csym
870 : INTEGER, DIMENSION(:, :, :), INTENT(OUT) :: srot
871 : INTEGER, INTENT(OUT) :: nrot
872 :
873 : INTEGER :: iop, jop, pass
874 : INTEGER, DIMENSION(3, 3) :: eye, frot, irot
875 : LOGICAL :: duplicate, identity
876 :
877 20 : CPASSERT(csym%symlib)
878 :
879 30388 : srot = 0
880 20 : eye = 0
881 20 : eye(1, 1) = 1
882 20 : eye(2, 2) = 1
883 20 : eye(3, 3) = 1
884 :
885 20 : nrot = 0
886 : ! Keep the identity first, matching the representative k-point operation.
887 60 : DO pass = 1, 2
888 4732 : DO iop = 1, csym%n_operations
889 60736 : irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
890 60736 : frot(1:3, 1:3) = TRANSPOSE(irot(1:3, 1:3))
891 10096 : identity = ALL(frot == eye)
892 4672 : IF (pass == 1 .AND. .NOT. identity) CYCLE
893 2392 : IF (pass == 2 .AND. identity) CYCLE
894 :
895 2336 : duplicate = .FALSE.
896 56528 : DO jop = 1, nrot
897 140872 : IF (ALL(frot == srot(:, :, jop))) THEN
898 : duplicate = .TRUE.
899 : EXIT
900 : END IF
901 : END DO
902 2336 : IF (duplicate) CYCLE
903 :
904 608 : nrot = nrot + 1
905 9728 : srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
906 : END DO
907 : END DO
908 :
909 20 : IF (nrot == 0) CALL cp_abort(__LOCATION__, "SPGLIB did not return usable symmetry rotations")
910 :
911 20 : END SUBROUTINE setup_spglib_reduction_rotations
912 :
913 : ! **************************************************************************************************
914 : !> \brief Determine the atom permutation generated by a SPGLIB space-group operation
915 : !> \param csym ...
916 : !> \param rot integer rotation in fractional coordinates
917 : !> \param trans fractional translation
918 : !> \param perm atom permutation
919 : !> \param valid whether all atoms were mapped
920 : ! **************************************************************************************************
921 16776 : SUBROUTINE spglib_atom_permutation(csym, rot, trans, perm, valid)
922 : TYPE(csym_type) :: csym
923 : INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
924 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: trans
925 : INTEGER, DIMENSION(:), INTENT(OUT) :: perm
926 : LOGICAL, INTENT(OUT) :: valid
927 :
928 : INTEGER :: i, j, nat
929 : LOGICAL :: found
930 16776 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
931 : REAL(KIND=dp) :: eps
932 : REAL(KIND=dp), DIMENSION(3) :: diff, spos
933 : REAL(KIND=dp), DIMENSION(3, 3) :: rfrac
934 :
935 16776 : nat = csym%nat
936 16776 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
937 218088 : rfrac(1:3, 1:3) = REAL(rot(1:3, 1:3), KIND=dp)
938 50328 : ALLOCATE (used(nat))
939 16776 : used = .FALSE.
940 119088 : perm = 0
941 16776 : valid = .TRUE.
942 :
943 119088 : DO i = 1, nat
944 1636992 : spos(1:3) = MATMUL(rfrac(1:3, 1:3), csym%scoord(1:3, i)) + trans(1:3)
945 459292 : found = .FALSE.
946 459292 : DO j = 1, nat
947 459292 : IF (used(j)) CYCLE
948 251210 : IF (csym%atype(i) /= csym%atype(j)) CYCLE
949 1004840 : diff(1:3) = spos(1:3) - csym%scoord(1:3, j)
950 1004840 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
951 577170 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
952 102312 : perm(i) = j
953 102312 : used(j) = .TRUE.
954 : found = .TRUE.
955 : EXIT
956 : END IF
957 : END DO
958 16776 : IF (.NOT. found) THEN
959 0 : valid = .FALSE.
960 0 : EXIT
961 : END IF
962 : END DO
963 :
964 16776 : DEALLOCATE (used)
965 :
966 16776 : END SUBROUTINE spglib_atom_permutation
967 :
968 : ! **************************************************************************************************
969 : !> \brief Reduce a k-point mesh with SPGLIB direct-space operations
970 : !> \param csym ...
971 : !> \param xkp full k-point mesh in reciprocal lattice coordinates
972 : !> \param wkp reduced k-point weights
973 : !> \param kpop symmetry operation mapping the representative k-point to a mesh point
974 : !> \param srot integer rotations in fractional coordinates
975 : !> \param nrot number of stored rotations
976 : ! **************************************************************************************************
977 994 : SUBROUTINE reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
978 : TYPE(csym_type) :: csym
979 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
980 : REAL(KIND=dp), DIMENSION(:) :: wkp
981 : INTEGER, DIMENSION(:) :: kpop
982 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
983 : INTEGER, INTENT(IN) :: nrot
984 :
985 : INTEGER :: i, iop, isign, j, kr, nkpts
986 994 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kscore, op_score
987 994 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: krot
988 : REAL(KIND=dp) :: eps
989 : REAL(KIND=dp), DIMENSION(3) :: rr
990 :
991 994 : nkpts = SIZE(wkp)
992 994 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
993 6958 : ALLOCATE (kscore(nkpts), op_score(nrot), krot(3, 3, nrot))
994 :
995 16906 : DO iop = 1, nrot
996 15912 : kr = csym%ibrot(iop)
997 206856 : krot(1:3, 1:3, iop) = reciprocal_rotation(srot(1:3, 1:3, kr))
998 16906 : op_score(iop) = spglib_operation_score(csym, iop, srot(1:3, 1:3, kr))
999 : END DO
1000 :
1001 14048 : wkp = 0.0_dp
1002 14048 : kpop = 0
1003 14048 : csym%kplink(1, :) = 0
1004 14048 : kscore = HUGE(0)
1005 :
1006 14048 : DO i = 1, nkpts
1007 13054 : IF (csym%kplink(1, i) /= 0) CYCLE
1008 :
1009 2436 : csym%kplink(1, i) = i
1010 2436 : wkp(i) = 1.0_dp
1011 2436 : kpop(i) = 1
1012 2436 : kscore(i) = 0
1013 :
1014 43584 : DO iop = 1, nrot
1015 133516 : DO isign = 1, 2
1016 2007700 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3, iop), KIND=dp), xkp(1:3, i))
1017 80308 : kr = csym%ibrot(iop)
1018 80308 : IF (isign == 2) THEN
1019 160616 : rr(1:3) = -rr(1:3)
1020 40154 : kr = -kr
1021 : END IF
1022 :
1023 80308 : j = regular_kpoint_index(rr, xkp(1:3, 1), csym%mesh, eps)
1024 80308 : IF (j == 0) CYCLE
1025 113566 : IF (csym%kplink(1, j) == 0) THEN
1026 10618 : csym%kplink(1, j) = i
1027 10618 : wkp(i) = wkp(i) + 1.0_dp
1028 10618 : kpop(j) = kr
1029 10618 : kscore(j) = op_score(iop)
1030 : ELSE
1031 62794 : CPASSERT(csym%kplink(1, j) == i)
1032 62794 : IF (op_score(iop) < kscore(j)) THEN
1033 3516 : kpop(j) = kr
1034 3516 : kscore(j) = op_score(iop)
1035 : END IF
1036 : END IF
1037 : END DO
1038 : END DO
1039 : END DO
1040 :
1041 14048 : DO i = 1, nkpts
1042 13054 : CPASSERT(csym%kplink(1, i) /= 0)
1043 14048 : CPASSERT(kpop(i) /= 0)
1044 : END DO
1045 994 : DEALLOCATE (kscore, op_score, krot)
1046 :
1047 994 : END SUBROUTINE reduce_spglib_kpoint_mesh
1048 :
1049 : ! **************************************************************************************************
1050 : !> \brief Reduce an explicit k-point set by inversion/time-reversal.
1051 : !> \param csym ...
1052 : !> \param xkp_full explicit k-point coordinates
1053 : !> \param wkp_full explicit k-point weights
1054 : ! **************************************************************************************************
1055 0 : SUBROUTINE reduce_general_inversion(csym, xkp_full, wkp_full)
1056 : TYPE(csym_type) :: csym
1057 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1058 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_full
1059 :
1060 : INTEGER :: i, j, nfull, nred
1061 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1062 0 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1063 : REAL(KIND=dp) :: eps
1064 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1065 : REAL(KIND=dp), DIMENSION(3) :: diff
1066 :
1067 0 : nfull = SIZE(wkp_full)
1068 0 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1069 0 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1070 0 : used = .FALSE.
1071 0 : rep = 0
1072 0 : wred = 0.0_dp
1073 0 : nred = 0
1074 :
1075 0 : DO i = 1, nfull
1076 0 : IF (used(i)) CYCLE
1077 0 : nred = nred + 1
1078 0 : rep(nred) = i
1079 0 : used(i) = .TRUE.
1080 0 : csym%kplink(1, i) = i
1081 0 : csym%kpop(i) = 1
1082 0 : wred(nred) = wkp_full(i)
1083 0 : DO j = i + 1, nfull
1084 0 : IF (used(j)) CYCLE
1085 0 : diff(1:3) = xkp_full(1:3, j) + xkp_full(1:3, i)
1086 0 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1087 0 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1088 0 : IF (ABS(wkp_full(j) - wkp_full(i)) > eps) THEN
1089 : CALL cp_abort(__LOCATION__, &
1090 : "KPOINTS%INVERSION_SYMMETRY_ONLY with SCHEME GENERAL requires "// &
1091 0 : "equal weights for inversion-related k-points.")
1092 : END IF
1093 0 : used(j) = .TRUE.
1094 0 : csym%kplink(1, j) = i
1095 0 : csym%kpop(j) = -1
1096 0 : wred(nred) = wred(nred) + wkp_full(j)
1097 : END IF
1098 : END DO
1099 : END DO
1100 :
1101 0 : csym%nkpoint = nred
1102 0 : ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1103 0 : DO i = 1, nred
1104 0 : csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1105 0 : csym%wkpoint(i) = wred(i)
1106 : END DO
1107 0 : DO i = 1, nfull
1108 0 : DO j = 1, nred
1109 0 : IF (csym%kplink(1, i) == rep(j)) THEN
1110 0 : csym%kplink(2, i) = j
1111 0 : EXIT
1112 : END IF
1113 : END DO
1114 : END DO
1115 :
1116 0 : DEALLOCATE (rep, used, wred)
1117 :
1118 0 : END SUBROUTINE reduce_general_inversion
1119 :
1120 : ! **************************************************************************************************
1121 : !> \brief Reduce an explicit k-point set with K290 symmetry operations.
1122 : !> \param csym ...
1123 : !> \param xkp_full explicit k-point coordinates
1124 : ! **************************************************************************************************
1125 4 : SUBROUTINE reduce_general_k290(csym, xkp_full)
1126 : TYPE(csym_type) :: csym
1127 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1128 :
1129 : INTEGER :: i, ibsign, iop, j, kr, nfull, nred
1130 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1131 : LOGICAL :: found
1132 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1133 : REAL(KIND=dp) :: alat, eps
1134 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1135 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, diff, rr, wcart
1136 :
1137 4 : nfull = SIZE(xkp_full, 2)
1138 4 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1139 4 : CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1140 :
1141 24 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1142 4 : used = .FALSE.
1143 4 : rep = 0
1144 4 : wred = 0.0_dp
1145 4 : nred = 0
1146 :
1147 36 : DO i = 1, nfull
1148 32 : IF (used(i)) CYCLE
1149 4 : nred = nred + 1
1150 4 : rep(nred) = i
1151 4 : used(i) = .TRUE.
1152 4 : csym%kplink(1, i) = i
1153 4 : csym%kpop(i) = 1
1154 4 : wred(nred) = 1.0_dp
1155 :
1156 200 : DO iop = 1, csym%nrtot
1157 608 : DO ibsign = 1, 2
1158 384 : kr = csym%ibrot(iop)
1159 : wcart(1:3) = alat*(xkp_full(1, i)*b1(1:3) + &
1160 : xkp_full(2, i)*b2(1:3) + &
1161 1536 : xkp_full(3, i)*b3(1:3))
1162 1536 : wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1163 384 : IF (ibsign == 2) THEN
1164 768 : wcart(1:3) = -wcart(1:3)
1165 192 : kr = -kr
1166 : END IF
1167 1536 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1168 1536 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1169 1536 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1170 :
1171 384 : found = .FALSE.
1172 1728 : DO j = 1, nfull
1173 6912 : diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1174 6912 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1175 3648 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1176 384 : found = .TRUE.
1177 384 : IF (.NOT. used(j)) THEN
1178 28 : used(j) = .TRUE.
1179 28 : csym%kplink(1, j) = i
1180 28 : csym%kpop(j) = kr
1181 28 : wred(nred) = wred(nred) + 1.0_dp
1182 : ELSE
1183 356 : CPASSERT(csym%kplink(1, j) == i)
1184 : END IF
1185 : EXIT
1186 : END IF
1187 : END DO
1188 192 : IF (.NOT. found) THEN
1189 : CALL cp_abort(__LOCATION__, &
1190 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1191 0 : "to be closed under the K290 symmetry operations.")
1192 : END IF
1193 : END DO
1194 : END DO
1195 : END DO
1196 :
1197 4 : CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1198 :
1199 4 : DEALLOCATE (rep, used, wred)
1200 :
1201 4 : END SUBROUTINE reduce_general_k290
1202 :
1203 : ! **************************************************************************************************
1204 : !> \brief Reduce an explicit k-point set with SPGLIB symmetry operations.
1205 : !> \param csym ...
1206 : !> \param xkp_full explicit k-point coordinates
1207 : ! **************************************************************************************************
1208 18 : SUBROUTINE reduce_general_spglib(csym, xkp_full)
1209 : TYPE(csym_type) :: csym
1210 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1211 :
1212 : INTEGER :: i, iop, isign, j, kr, nfull, nred, nrot
1213 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1214 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
1215 : INTEGER, DIMENSION(3, 3) :: krot
1216 : LOGICAL :: found
1217 18 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1218 : REAL(KIND=dp) :: eps
1219 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1220 : REAL(KIND=dp), DIMENSION(3) :: diff, rr
1221 :
1222 18 : nfull = SIZE(xkp_full, 2)
1223 18 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1224 54 : ALLOCATE (srot(3, 3, csym%n_operations))
1225 18 : CALL setup_spglib_operations(csym, srot, nrot)
1226 :
1227 108 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1228 18 : used = .FALSE.
1229 18 : rep = 0
1230 18 : wred = 0.0_dp
1231 18 : nred = 0
1232 :
1233 162 : DO i = 1, nfull
1234 144 : IF (used(i)) CYCLE
1235 18 : nred = nred + 1
1236 18 : rep(nred) = i
1237 18 : used(i) = .TRUE.
1238 18 : csym%kplink(1, i) = i
1239 18 : csym%kpop(i) = 1
1240 18 : wred(nred) = 1.0_dp
1241 :
1242 900 : DO iop = 1, nrot
1243 864 : kr = csym%ibrot(iop)
1244 864 : krot = reciprocal_rotation(srot(:, :, kr))
1245 2736 : DO isign = 1, 2
1246 43200 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp_full(1:3, i))
1247 1728 : IF (isign == 2) THEN
1248 3456 : rr(1:3) = -rr(1:3)
1249 864 : kr = -csym%ibrot(iop)
1250 : ELSE
1251 864 : kr = csym%ibrot(iop)
1252 : END IF
1253 :
1254 1728 : found = .FALSE.
1255 7776 : DO j = 1, nfull
1256 31104 : diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1257 31104 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1258 16416 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1259 1728 : found = .TRUE.
1260 1728 : IF (.NOT. used(j)) THEN
1261 126 : used(j) = .TRUE.
1262 126 : csym%kplink(1, j) = i
1263 126 : csym%kpop(j) = kr
1264 126 : wred(nred) = wred(nred) + 1.0_dp
1265 : ELSE
1266 1602 : CPASSERT(csym%kplink(1, j) == i)
1267 : END IF
1268 : EXIT
1269 : END IF
1270 : END DO
1271 864 : IF (.NOT. found) THEN
1272 : CALL cp_abort(__LOCATION__, &
1273 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1274 0 : "to be closed under the requested symmetry operations.")
1275 : END IF
1276 : END DO
1277 : END DO
1278 : END DO
1279 :
1280 18 : CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1281 :
1282 18 : DEALLOCATE (rep, srot, used, wred)
1283 :
1284 18 : END SUBROUTINE reduce_general_spglib
1285 :
1286 : ! **************************************************************************************************
1287 : !> \brief Reduce an explicit k-point set with SPGLIB rotations and K290 operations.
1288 : !> \param csym ...
1289 : !> \param xkp_full explicit k-point coordinates
1290 : ! **************************************************************************************************
1291 4 : SUBROUTINE reduce_general_spglib_k290(csym, xkp_full)
1292 : TYPE(csym_type) :: csym
1293 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1294 :
1295 : INTEGER :: i, iop, isign, j, k290_op, nfull, nred, &
1296 : nrot, nskipped
1297 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1298 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
1299 : INTEGER, DIMENSION(3, 3) :: krot
1300 : LOGICAL :: found, valid
1301 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1302 : REAL(KIND=dp) :: alat, eps
1303 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1304 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, diff, rr
1305 :
1306 4 : nfull = SIZE(xkp_full, 2)
1307 4 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1308 4 : CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1309 12 : ALLOCATE (srot(3, 3, csym%n_operations))
1310 4 : CALL setup_spglib_reduction_rotations(csym, srot, nrot)
1311 :
1312 24 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1313 4 : used = .FALSE.
1314 4 : rep = 0
1315 4 : wred = 0.0_dp
1316 4 : nred = 0
1317 4 : nskipped = 0
1318 :
1319 36 : DO i = 1, nfull
1320 32 : IF (used(i)) CYCLE
1321 4 : nred = nred + 1
1322 4 : rep(nred) = i
1323 4 : used(i) = .TRUE.
1324 4 : csym%kplink(1, i) = i
1325 4 : csym%kpop(i) = 1
1326 4 : wred(nred) = 1.0_dp
1327 :
1328 200 : DO iop = 1, nrot
1329 192 : krot = reciprocal_rotation(srot(:, :, iop))
1330 608 : DO isign = 1, 2
1331 9600 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp_full(1:3, i))
1332 960 : IF (isign == 2) rr(1:3) = -rr(1:3)
1333 :
1334 1728 : found = .FALSE.
1335 1728 : DO j = 1, nfull
1336 6912 : diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1337 6912 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1338 3648 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1339 384 : found = .TRUE.
1340 : CALL find_k290_kpoint_operation(csym, xkp_full(1:3, i), xkp_full(1:3, j), &
1341 : a1, a2, a3, b1, b2, b3, alat, &
1342 384 : k290_op, valid)
1343 384 : IF (.NOT. valid) THEN
1344 0 : nskipped = nskipped + 1
1345 : EXIT
1346 : END IF
1347 384 : IF (.NOT. used(j)) THEN
1348 28 : used(j) = .TRUE.
1349 28 : csym%kplink(1, j) = i
1350 28 : csym%kpop(j) = k290_op
1351 28 : wred(nred) = wred(nred) + 1.0_dp
1352 : ELSE
1353 356 : CPASSERT(csym%kplink(1, j) == i)
1354 : END IF
1355 : EXIT
1356 : END IF
1357 : END DO
1358 192 : IF (.NOT. found) THEN
1359 : CALL cp_abort(__LOCATION__, &
1360 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1361 0 : "to be closed under the SPGLIB symmetry operations.")
1362 : END IF
1363 : END DO
1364 : END DO
1365 : END DO
1366 :
1367 4 : IF (nskipped > 0) THEN
1368 : CALL cp_warn(__LOCATION__, &
1369 : "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1370 0 : "the GENERAL k-point set was reduced only by the compatible mappings.")
1371 : END IF
1372 :
1373 4 : CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1374 :
1375 4 : DEALLOCATE (rep, srot, used, wred)
1376 :
1377 8 : END SUBROUTINE reduce_general_spglib_k290
1378 :
1379 : ! **************************************************************************************************
1380 : !> \brief Store reduced GENERAL k-point representatives.
1381 : !> \param csym ...
1382 : !> \param xkp_full explicit k-point coordinates
1383 : !> \param rep representative indices
1384 : !> \param wred representative multiplicities
1385 : !> \param nred number of reduced representatives
1386 : ! **************************************************************************************************
1387 26 : SUBROUTINE store_general_reduction(csym, xkp_full, rep, wred, nred)
1388 : TYPE(csym_type) :: csym
1389 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1390 : INTEGER, DIMENSION(:), INTENT(IN) :: rep
1391 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wred
1392 : INTEGER, INTENT(IN) :: nred
1393 :
1394 : INTEGER :: i, j, nfull
1395 :
1396 26 : nfull = SIZE(xkp_full, 2)
1397 26 : csym%nkpoint = nred
1398 130 : ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1399 52 : DO i = 1, nred
1400 104 : csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1401 52 : csym%wkpoint(i) = wred(i)
1402 : END DO
1403 234 : DO i = 1, nfull
1404 208 : DO j = 1, nred
1405 208 : IF (csym%kplink(1, i) == rep(j)) THEN
1406 208 : csym%kplink(2, i) = j
1407 208 : EXIT
1408 : END IF
1409 : END DO
1410 234 : CPASSERT(csym%kplink(2, i) /= 0)
1411 : END DO
1412 :
1413 26 : END SUBROUTINE store_general_reduction
1414 :
1415 : ! **************************************************************************************************
1416 : !> \brief Reduce a k-point mesh with SPGLIB rotations and K290 operations
1417 : !> \param csym ...
1418 : !> \param xkp full k-point mesh in reciprocal lattice coordinates
1419 : !> \param wkp reduced k-point weights
1420 : !> \param kpop K290 operation mapping the representative k-point to a mesh point
1421 : !> \param srot SPGLIB integer rotations in fractional coordinates
1422 : !> \param nrot number of stored rotations
1423 : !> \param a1 first lattice vector
1424 : !> \param a2 second lattice vector
1425 : !> \param a3 third lattice vector
1426 : !> \param b1 first reciprocal lattice vector
1427 : !> \param b2 second reciprocal lattice vector
1428 : !> \param b3 third reciprocal lattice vector
1429 : !> \param alat lattice scaling used by K290
1430 : ! **************************************************************************************************
1431 16 : SUBROUTINE reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
1432 : a1, a2, a3, b1, b2, b3, alat)
1433 : TYPE(csym_type) :: csym
1434 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1435 : REAL(KIND=dp), DIMENSION(:) :: wkp
1436 : INTEGER, DIMENSION(:) :: kpop
1437 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
1438 : INTEGER, INTENT(IN) :: nrot
1439 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1440 : REAL(KIND=dp), INTENT(IN) :: alat
1441 :
1442 : INTEGER :: i, iop, isign, j, k290_op, nkpts, &
1443 : nskipped
1444 : INTEGER, DIMENSION(3, 3) :: krot
1445 : LOGICAL :: valid
1446 : REAL(KIND=dp) :: eps
1447 : REAL(KIND=dp), DIMENSION(3) :: rr
1448 :
1449 16 : nkpts = SIZE(wkp)
1450 16 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1451 16 : nskipped = 0
1452 :
1453 368 : wkp = 0.0_dp
1454 368 : kpop = 0
1455 368 : csym%kplink(1, :) = 0
1456 :
1457 368 : DO i = 1, nkpts
1458 352 : IF (csym%kplink(1, i) /= 0) CYCLE
1459 :
1460 64 : csym%kplink(1, i) = i
1461 64 : wkp(i) = 1.0_dp
1462 64 : kpop(i) = 1
1463 :
1464 688 : DO iop = 1, nrot
1465 608 : krot = reciprocal_rotation(srot(:, :, iop))
1466 2176 : DO isign = 1, 2
1467 30400 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp(1:3, i))
1468 3040 : IF (isign == 2) rr(1:3) = -rr(1:3)
1469 :
1470 1216 : j = regular_kpoint_index(rr, xkp(1:3, 1), csym%mesh, eps)
1471 1216 : IF (j == 0 .OR. j == i) CYCLE
1472 1024 : IF (csym%kplink(1, j) /= 0) THEN
1473 736 : CPASSERT(csym%kplink(1, j) == i)
1474 : CYCLE
1475 : END IF
1476 :
1477 : CALL find_k290_kpoint_operation(csym, xkp(1:3, i), xkp(1:3, j), &
1478 : a1, a2, a3, b1, b2, b3, alat, &
1479 288 : k290_op, valid)
1480 288 : IF (.NOT. valid) THEN
1481 0 : nskipped = nskipped + 1
1482 0 : CYCLE
1483 : END IF
1484 288 : csym%kplink(1, j) = i
1485 288 : wkp(i) = wkp(i) + 1.0_dp
1486 896 : kpop(j) = k290_op
1487 : END DO
1488 : END DO
1489 : END DO
1490 :
1491 368 : DO i = 1, nkpts
1492 352 : CPASSERT(csym%kplink(1, i) /= 0)
1493 368 : CPASSERT(kpop(i) /= 0)
1494 : END DO
1495 16 : IF (nskipped > 0) THEN
1496 : CALL cp_warn(__LOCATION__, &
1497 : "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1498 0 : "the mesh was reduced only by the compatible mappings.")
1499 : END IF
1500 :
1501 16 : END SUBROUTINE reduce_spglib_kpoint_mesh_k290
1502 :
1503 : ! **************************************************************************************************
1504 : !> \brief Find a K290 operation that maps one fractional k-point to another
1505 : !> \param csym ...
1506 : !> \param xref representative k-point
1507 : !> \param xtarget target k-point
1508 : !> \param a1 first lattice vector
1509 : !> \param a2 second lattice vector
1510 : !> \param a3 third lattice vector
1511 : !> \param b1 first reciprocal lattice vector
1512 : !> \param b2 second reciprocal lattice vector
1513 : !> \param b3 third reciprocal lattice vector
1514 : !> \param alat lattice scaling used by K290
1515 : !> \param k290_op K290 operation identifier
1516 : !> \param valid whether a matching K290 operation was found
1517 : ! **************************************************************************************************
1518 672 : SUBROUTINE find_k290_kpoint_operation(csym, xref, xtarget, a1, a2, a3, b1, b2, b3, alat, &
1519 : k290_op, valid)
1520 : TYPE(csym_type) :: csym
1521 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xref, xtarget, a1, a2, a3, b1, b2, b3
1522 : REAL(KIND=dp), INTENT(IN) :: alat
1523 : INTEGER, INTENT(OUT) :: k290_op
1524 : LOGICAL, INTENT(OUT) :: valid
1525 :
1526 : INTEGER :: ibsign, iop, kr
1527 : REAL(KIND=dp) :: eps
1528 : REAL(KIND=dp), DIMENSION(3) :: diff, rr, wcart
1529 :
1530 672 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1531 672 : k290_op = 0
1532 672 : valid = .FALSE.
1533 :
1534 1616 : DO iop = 1, csym%nrtot
1535 1616 : IF (iop > SIZE(csym%rt, 3)) CYCLE
1536 1616 : IF (csym%ibrot(iop) == 0) CYCLE
1537 3872 : DO ibsign = 1, 2
1538 11712 : wcart(1:3) = alat*(xref(1)*b1(1:3) + xref(2)*b2(1:3) + xref(3)*b3(1:3))
1539 11712 : wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1540 2928 : IF (ibsign == 2) THEN
1541 5248 : wcart(1:3) = -wcart(1:3)
1542 1312 : kr = -csym%ibrot(iop)
1543 : ELSE
1544 1616 : kr = csym%ibrot(iop)
1545 : END IF
1546 11712 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1547 11712 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1548 11712 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1549 :
1550 11712 : diff(1:3) = xtarget(1:3) - rr(1:3)
1551 11712 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1552 7056 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1553 672 : k290_op = kr
1554 672 : valid = .TRUE.
1555 672 : RETURN
1556 : END IF
1557 : END DO
1558 : END DO
1559 :
1560 : END SUBROUTINE find_k290_kpoint_operation
1561 :
1562 : ! **************************************************************************************************
1563 : !> \brief Score SPGLIB operations to choose stable atom transformations
1564 : !> \param csym ...
1565 : !> \param iop operation index
1566 : !> \param srot integer rotation in fractional coordinates
1567 : !> \return score, lower values are preferred
1568 : ! **************************************************************************************************
1569 15912 : FUNCTION spglib_operation_score(csym, iop, srot) RESULT(score)
1570 : TYPE(csym_type), INTENT(IN) :: csym
1571 : INTEGER, INTENT(IN) :: iop
1572 : INTEGER, DIMENSION(3, 3), INTENT(IN) :: srot
1573 : INTEGER :: score
1574 :
1575 : INTEGER :: i, nat
1576 : INTEGER, DIMENSION(3, 3) :: eye, r2
1577 : REAL(KIND=dp) :: eps
1578 :
1579 15912 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1580 15912 : nat = SIZE(csym%f0, 1)
1581 15912 : score = 0
1582 111312 : DO i = 1, nat
1583 111312 : IF (csym%f0(i, iop) /= i) score = score + 100
1584 : END DO
1585 45682 : IF (ANY(ABS(csym%vt(1:3, iop) - ANINT(csym%vt(1:3, iop))) > eps)) score = score + 10
1586 :
1587 15912 : eye = 0
1588 15912 : eye(1, 1) = 1
1589 15912 : eye(2, 2) = 1
1590 15912 : eye(3, 3) = 1
1591 636480 : r2(1:3, 1:3) = MATMUL(srot(1:3, 1:3), srot(1:3, 1:3))
1592 112720 : IF (ANY(r2(1:3, 1:3) /= eye(1:3, 1:3))) score = score + 1
1593 :
1594 15912 : END FUNCTION spglib_operation_score
1595 :
1596 : ! **************************************************************************************************
1597 : !> \brief Reciprocal-space rotation corresponding to a fractional direct-space rotation
1598 : !> \param rot direct-space rotation
1599 : !> \return reciprocal-space rotation
1600 : ! **************************************************************************************************
1601 17576 : FUNCTION reciprocal_rotation(rot) RESULT(krot)
1602 : INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
1603 : INTEGER, DIMENSION(3, 3) :: krot
1604 :
1605 : REAL(KIND=dp), DIMENSION(3, 3) :: rinv
1606 :
1607 228488 : rinv = inv_3x3(REAL(rot(1:3, 1:3), KIND=dp))
1608 228488 : krot(1:3, 1:3) = NINT(TRANSPOSE(rinv(1:3, 1:3)))
1609 :
1610 17576 : END FUNCTION reciprocal_rotation
1611 :
1612 : ! **************************************************************************************************
1613 : !> \brief Reduce a CP2K Monkhorst-Pack mesh using K290 symmetry operations
1614 : !> \param csym ...
1615 : !> \param xkp full k-point mesh in reciprocal lattice coordinates
1616 : !> \param wkp reduced k-point weights
1617 : !> \param kpop symmetry operation mapping the representative k-point to a mesh point
1618 : !> \param nc number of point group operations
1619 : !> \param ib K290 operation identifiers
1620 : !> \param r K290 rotation matrices
1621 : !> \param a1 first lattice vector
1622 : !> \param a2 second lattice vector
1623 : !> \param a3 third lattice vector
1624 : !> \param b1 first reciprocal lattice vector
1625 : !> \param b2 second reciprocal lattice vector
1626 : !> \param b3 third reciprocal lattice vector
1627 : !> \param alat lattice scaling used by K290
1628 : ! **************************************************************************************************
1629 700 : SUBROUTINE reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
1630 : TYPE(csym_type) :: csym
1631 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1632 : REAL(KIND=dp), DIMENSION(:) :: wkp
1633 : INTEGER, DIMENSION(:) :: kpop
1634 : INTEGER, INTENT(IN) :: nc
1635 : INTEGER, DIMENSION(48), INTENT(IN) :: ib
1636 : REAL(KIND=dp), DIMENSION(3, 3, 48), INTENT(IN) :: r
1637 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1638 : REAL(KIND=dp), INTENT(IN) :: alat
1639 :
1640 : INTEGER :: i, ibsign, iop, j, kr, nkpts, &
1641 : nskipped_overlaps
1642 : REAL(KIND=dp) :: eps
1643 : REAL(KIND=dp), DIMENSION(3) :: rr, wcart
1644 :
1645 700 : nkpts = SIZE(wkp)
1646 700 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1647 700 : nskipped_overlaps = 0
1648 :
1649 7146 : wkp = 0.0_dp
1650 7146 : kpop = 0
1651 7146 : csym%kplink(1, :) = 0
1652 :
1653 7146 : DO i = 1, nkpts
1654 6446 : IF (csym%kplink(1, i) /= 0) CYCLE
1655 :
1656 1684 : csym%kplink(1, i) = i
1657 1684 : wkp(i) = 1.0_dp
1658 1684 : kpop(i) = 1
1659 :
1660 25580 : DO iop = 1, nc
1661 76034 : DO ibsign = 1, 2
1662 46392 : kr = ib(iop)
1663 185568 : wcart(1:3) = alat*(xkp(1, i)*b1(1:3) + xkp(2, i)*b2(1:3) + xkp(3, i)*b3(1:3))
1664 185568 : wcart(1:3) = kp_apply_operation(wcart(1:3), r(1:3, 1:3, kr))
1665 46392 : IF (ibsign == 2) THEN
1666 92784 : wcart(1:3) = -wcart(1:3)
1667 23196 : kr = -kr
1668 : END IF
1669 185568 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1670 185568 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1671 185568 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1672 :
1673 46392 : j = regular_kpoint_index(rr, xkp(1:3, 1), csym%mesh, eps)
1674 : ! Some point-group operations are incompatible with the requested Monkhorst-Pack mesh.
1675 46392 : IF (j == 0) CYCLE
1676 68856 : IF (csym%kplink(1, j) == 0) THEN
1677 4762 : csym%kplink(1, j) = i
1678 4762 : wkp(i) = wkp(i) + 1.0_dp
1679 4762 : kpop(j) = kr
1680 40898 : ELSE IF (csym%kplink(1, j) /= i) THEN
1681 : ! Approximate K290 operation sets need not be closed for structures whose
1682 : ! coordinates lie close to several symmetry tolerances. Keep the existing
1683 : ! disjoint orbit instead of aborting or double-counting this mesh point.
1684 16 : nskipped_overlaps = nskipped_overlaps + 1
1685 : END IF
1686 : END DO
1687 : END DO
1688 : END DO
1689 :
1690 700 : IF (nskipped_overlaps > 0) THEN
1691 56 : wkp = 1.0_dp
1692 2 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
1693 56 : DO i = 1, nkpts
1694 56 : IF (wkp(i) > 0.0_dp) THEN
1695 28 : kpop(i) = 1
1696 : ELSE
1697 26 : kpop(i) = 2
1698 : END IF
1699 : END DO
1700 2 : csym%nrtot = 0
1701 2 : csym%inversion_only = .TRUE.
1702 : CALL cp_warn(__LOCATION__, &
1703 : "The K290 k-point operations produced overlapping, non-closed mesh orbits; "// &
1704 : "falling back to inversion/time-reversal symmetry. Use SYMMETRY_BACKEND SPGLIB "// &
1705 2 : "for full reduction by a closed crystallographic symmetry group.")
1706 : ELSE
1707 7090 : DO i = 1, nkpts
1708 6392 : CPASSERT(csym%kplink(1, i) /= 0)
1709 7090 : CPASSERT(kpop(i) /= 0)
1710 : END DO
1711 : END IF
1712 :
1713 700 : END SUBROUTINE reduce_kpoint_mesh
1714 :
1715 : ! **************************************************************************************************
1716 : !> \brief Return the flat index of a point on a regular Monkhorst-Pack mesh.
1717 : !> \param x target k-point in reciprocal lattice coordinates
1718 : !> \param origin first point of the generated mesh
1719 : !> \param mesh mesh dimensions
1720 : !> \param eps coordinate tolerance
1721 : !> \return flat mesh index, or zero if x is not on the mesh
1722 : ! **************************************************************************************************
1723 127916 : PURE FUNCTION regular_kpoint_index(x, origin, mesh, eps) RESULT(index)
1724 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: x, origin
1725 : INTEGER, DIMENSION(3), INTENT(IN) :: mesh
1726 : REAL(KIND=dp), INTENT(IN) :: eps
1727 : INTEGER :: index
1728 :
1729 : INTEGER :: idim
1730 : INTEGER, DIMENSION(3) :: icoord, istep
1731 : REAL(KIND=dp) :: scaled
1732 :
1733 127916 : index = 0
1734 511664 : IF (ANY(mesh <= 0)) RETURN
1735 494764 : DO idim = 1, 3
1736 374476 : scaled = (x(idim) - origin(idim))*REAL(mesh(idim), KIND=dp)
1737 374476 : istep(idim) = NINT(scaled)
1738 374476 : IF (ABS(scaled - REAL(istep(idim), KIND=dp)) > eps*REAL(mesh(idim), KIND=dp)) RETURN
1739 487136 : icoord(idim) = MODULO(istep(idim), mesh(idim)) + 1
1740 : END DO
1741 :
1742 120288 : index = ((icoord(1) - 1)*mesh(2) + icoord(2) - 1)*mesh(3) + icoord(3)
1743 :
1744 120288 : END FUNCTION regular_kpoint_index
1745 :
1746 : ! **************************************************************************************************
1747 : !> \brief ...
1748 : !> \param nk ...
1749 : !> \param xkp ...
1750 : !> \param wkp ...
1751 : !> \param shift ...
1752 : !> \param gamma_centered ...
1753 : ! **************************************************************************************************
1754 3430 : SUBROUTINE full_grid_gen(nk, xkp, wkp, shift, gamma_centered)
1755 : INTEGER, INTENT(IN) :: nk(3)
1756 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
1757 : REAL(KIND=dp), DIMENSION(:) :: wkp
1758 : REAL(KIND=dp), INTENT(IN) :: shift(3)
1759 : LOGICAL, INTENT(IN), OPTIONAL :: gamma_centered
1760 :
1761 : INTEGER :: i, idim, ix, iy, iz
1762 : INTEGER, DIMENSION(3) :: ik
1763 : LOGICAL :: gamma_mesh
1764 : REAL(KIND=dp) :: kpt_latt(3)
1765 :
1766 3430 : IF (PRESENT(gamma_centered)) THEN
1767 3430 : gamma_mesh = gamma_centered
1768 : ELSE
1769 : gamma_mesh = .FALSE.
1770 : END IF
1771 :
1772 37898 : wkp = 0.0_dp
1773 3430 : i = 0
1774 10212 : DO ix = 1, nk(1)
1775 24560 : DO iy = 1, nk(2)
1776 55598 : DO iz = 1, nk(3)
1777 34468 : i = i + 1
1778 34468 : ik(1) = ix
1779 34468 : ik(2) = iy
1780 34468 : ik(3) = iz
1781 137872 : DO idim = 1, 3
1782 137872 : IF (gamma_mesh .AND. MOD(nk(idim), 2) == 0) THEN
1783 : kpt_latt(idim) = REAL(2*ik(idim) - nk(idim), KIND=dp)/ &
1784 2148 : (2._dp*REAL(nk(idim), KIND=dp))
1785 : ELSE
1786 : kpt_latt(idim) = REAL(2*ik(idim) - nk(idim) - 1, KIND=dp)/ &
1787 101256 : (2._dp*REAL(nk(idim), KIND=dp))
1788 : END IF
1789 : END DO
1790 137872 : xkp(1:3, i) = kpt_latt(1:3)
1791 48816 : wkp(i) = 1.0_dp
1792 : END DO
1793 : END DO
1794 : END DO
1795 37898 : DO i = 1, nk(1)*nk(2)*nk(3)
1796 141302 : xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
1797 : END DO
1798 :
1799 3430 : END SUBROUTINE full_grid_gen
1800 :
1801 : ! **************************************************************************************************
1802 : !> \brief ...
1803 : !> \param xkp ...
1804 : !> \param wkp ...
1805 : !> \param link ...
1806 : ! **************************************************************************************************
1807 1620 : SUBROUTINE inversion_symm(xkp, wkp, link)
1808 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
1809 : REAL(KIND=dp), DIMENSION(:) :: wkp
1810 : INTEGER, DIMENSION(:) :: link
1811 :
1812 : INTEGER :: i, j, nkpts
1813 : REAL(KIND=dp), DIMENSION(3) :: diff
1814 :
1815 1620 : nkpts = SIZE(wkp, 1)
1816 :
1817 15244 : link(:) = 0
1818 15244 : DO i = 1, nkpts
1819 13624 : IF (link(i) == 0) link(i) = i
1820 171284 : DO j = i + 1, nkpts
1821 162564 : IF (wkp(j) == 0) CYCLE
1822 444112 : diff(1:3) = xkp(1:3, i) + xkp(1:3, j)
1823 444112 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1824 166770 : IF (ALL(ABS(diff(1:3)) < 1.e-12_dp)) THEN
1825 6524 : wkp(i) = wkp(i) + wkp(j)
1826 6524 : wkp(j) = 0.0_dp
1827 6524 : link(j) = i
1828 6524 : EXIT
1829 : END IF
1830 : END DO
1831 : END DO
1832 :
1833 1620 : END SUBROUTINE inversion_symm
1834 :
1835 : ! **************************************************************************************************
1836 : !> \brief ...
1837 : !> \param x ...
1838 : !> \param r ...
1839 : !> \return ...
1840 : ! **************************************************************************************************
1841 49704 : FUNCTION kp_apply_operation(x, r) RESULT(y)
1842 : REAL(KIND=dp), INTENT(IN) :: x(3), r(3, 3)
1843 : REAL(KIND=dp) :: y(3)
1844 :
1845 49704 : y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
1846 49704 : y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
1847 49704 : y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
1848 :
1849 49704 : END FUNCTION kp_apply_operation
1850 :
1851 : ! **************************************************************************************************
1852 : !> \brief ...
1853 : !> \param csym ...
1854 : ! **************************************************************************************************
1855 3468 : SUBROUTINE print_crys_symmetry(csym)
1856 : TYPE(csym_type) :: csym
1857 :
1858 : INTEGER :: i, iunit, j, plevel
1859 :
1860 3468 : iunit = csym%punit
1861 3468 : IF (iunit >= 0) THEN
1862 1469 : plevel = csym%plevel
1863 1469 : WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
1864 1469 : IF (csym%symlib) THEN
1865 1467 : WRITE (iunit, '(A,T71,A10)') " International Symbol: ", ADJUSTR(TRIM(csym%international_symbol))
1866 1467 : WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", ADJUSTR(TRIM(csym%pointgroup_symbol))
1867 1467 : WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", ADJUSTR(TRIM(csym%schoenflies))
1868 : !
1869 1467 : WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
1870 1467 : IF (plevel > 0) THEN
1871 0 : DO i = 1, csym%n_operations
1872 : WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
1873 0 : " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
1874 0 : WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
1875 : END DO
1876 : END IF
1877 : ELSE
1878 2 : IF (csym%spglib_requested) THEN
1879 1 : WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not available"
1880 : ELSE
1881 1 : WRITE (iunit, "(T2,A)") "SPGLIB Crystal Symmetry Information was not requested"
1882 : END IF
1883 : END IF
1884 : END IF
1885 :
1886 3468 : END SUBROUTINE print_crys_symmetry
1887 :
1888 : ! **************************************************************************************************
1889 : !> \brief ...
1890 : !> \param csym ...
1891 : ! **************************************************************************************************
1892 3178 : SUBROUTINE print_kp_symmetry(csym)
1893 : TYPE(csym_type), INTENT(IN) :: csym
1894 :
1895 : INTEGER :: i, iunit, nat, nmesh, plevel
1896 :
1897 3178 : iunit = csym%punit
1898 3178 : IF (iunit >= 0) THEN
1899 1179 : plevel = csym%plevel
1900 1179 : WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
1901 1179 : WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
1902 1179 : WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
1903 4569 : DO i = 1, csym%nkpoint
1904 4569 : WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), NINT(csym%wkpoint(i))
1905 : END DO
1906 1179 : nmesh = csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
1907 1179 : IF (nmesh > 0) THEN
1908 1171 : WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
1909 : ELSE
1910 8 : nmesh = SIZE(csym%kpmesh, 2)
1911 8 : WRITE (iunit, '(/,A,T70,I10)') " Explicit K-point Set: ", nmesh
1912 : END IF
1913 1179 : WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
1914 11262 : DO i = 1, nmesh
1915 10083 : WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
1916 21345 : csym%kplink(1:2, i), csym%kpop(i)
1917 : END DO
1918 1179 : IF (csym%nrtot > 0) THEN
1919 662 : WRITE (iunit, '(/,A)') " Atom Transformation Table"
1920 662 : nat = SIZE(csym%f0, 1)
1921 7568 : DO i = 1, csym%nrtot
1922 7568 : WRITE (iunit, '(T10,A,I5,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
1923 : END DO
1924 : END IF
1925 : END IF
1926 :
1927 3178 : END SUBROUTINE print_kp_symmetry
1928 :
1929 0 : END MODULE cryssym
|
