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 3808 : SUBROUTINE release_csym_type(csym)
89 : TYPE(csym_type) :: csym
90 :
91 3808 : IF (ALLOCATED(csym%rotations)) THEN
92 3497 : DEALLOCATE (csym%rotations)
93 : END IF
94 3808 : IF (ALLOCATED(csym%translations)) THEN
95 3497 : DEALLOCATE (csym%translations)
96 : END IF
97 3808 : IF (ALLOCATED(csym%atype)) THEN
98 3808 : DEALLOCATE (csym%atype)
99 : END IF
100 3808 : IF (ALLOCATED(csym%scoord)) THEN
101 3808 : DEALLOCATE (csym%scoord)
102 : END IF
103 3808 : IF (ALLOCATED(csym%xkpoint)) THEN
104 3518 : DEALLOCATE (csym%xkpoint)
105 : END IF
106 3808 : IF (ALLOCATED(csym%wkpoint)) THEN
107 3518 : DEALLOCATE (csym%wkpoint)
108 : END IF
109 3808 : IF (ALLOCATED(csym%kpmesh)) THEN
110 3518 : DEALLOCATE (csym%kpmesh)
111 : END IF
112 3808 : IF (ALLOCATED(csym%kplink)) THEN
113 3518 : DEALLOCATE (csym%kplink)
114 : END IF
115 3808 : IF (ALLOCATED(csym%kpop)) THEN
116 3518 : DEALLOCATE (csym%kpop)
117 : END IF
118 3808 : IF (ALLOCATED(csym%rt)) THEN
119 3518 : DEALLOCATE (csym%rt)
120 : END IF
121 3808 : IF (ALLOCATED(csym%vt)) THEN
122 3518 : DEALLOCATE (csym%vt)
123 : END IF
124 3808 : IF (ALLOCATED(csym%f0)) THEN
125 3518 : DEALLOCATE (csym%f0)
126 : END IF
127 3808 : IF (ALLOCATED(csym%ibrot)) THEN
128 3518 : DEALLOCATE (csym%ibrot)
129 : END IF
130 :
131 3808 : 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 3808 : 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 3808 : CALL timeset(routineN, handle)
159 :
160 : !..total number of atoms
161 3808 : nat = SIZE(scoor, 2)
162 3808 : csym%nat = nat
163 :
164 : ! output unit
165 3808 : IF (PRESENT(iounit)) THEN
166 3796 : csym%punit = iounit
167 : ELSE
168 12 : csym%punit = -1
169 : END IF
170 :
171 : ! accuracy for symmetry
172 3808 : IF (PRESENT(delta)) THEN
173 3808 : csym%delta = delta
174 : ELSE
175 0 : csym%delta = 1.e-6_dp
176 : END IF
177 :
178 : !..set cell values
179 49504 : csym%hmat = hmat
180 :
181 : ! atom types
182 11424 : ALLOCATE (csym%atype(nat))
183 25415 : csym%atype(1:nat) = types(1:nat)
184 :
185 : ! scaled coordinates
186 11424 : ALLOCATE (csym%scoord(3, nat))
187 90236 : csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
188 :
189 3808 : csym%n_operations = 0
190 :
191 : !..try spglib
192 3808 : my_use_spglib = .TRUE.
193 3808 : IF (PRESENT(use_spglib)) my_use_spglib = use_spglib
194 290 : csym%spglib_requested = my_use_spglib
195 3518 : IF (.NOT. my_use_spglib) THEN
196 : spglib = .FALSE.
197 : ELSE
198 3498 : major = spg_get_major_version()
199 3498 : minor = spg_get_minor_version()
200 3498 : micro = spg_get_micro_version()
201 3498 : IF (major == 0) THEN
202 0 : CALL cp_warn(__LOCATION__, "Symmetry library SPGLIB not available")
203 0 : spglib = .FALSE.
204 : ELSE
205 3498 : spglib = .TRUE.
206 3498 : CALL cite_reference(Togo2018)
207 3498 : ierr = spg_get_international(csym%international_symbol, TRANSPOSE(hmat), scoor, types, nat, delta)
208 3498 : IF (ierr == 0) THEN
209 1 : CALL cp_warn(__LOCATION__, "Symmetry Library SPGLIB failed")
210 1 : spglib = .FALSE.
211 : ELSE
212 3497 : nop = spg_get_multiplicity(TRANSPOSE(hmat), scoor, types, nat, delta)
213 17485 : ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
214 3497 : csym%n_operations = nop
215 : ierr = spg_get_symmetry(csym%rotations, csym%translations, nop, &
216 3497 : TRANSPOSE(hmat), scoor, types, nat, delta)
217 : ! Schoenflies Symbol
218 3497 : csym%schoenflies = ' '
219 3497 : ierr = spg_get_schoenflies(csym%schoenflies(1:7), TRANSPOSE(hmat), scoor, types, nat, delta)
220 : ! Point Group
221 3497 : csym%pointgroup_symbol = ' '
222 3497 : tra_mat = 0
223 : ierr = spg_get_pointgroup(csym%pointgroup_symbol, tra_mat, &
224 3497 : csym%rotations, csym%n_operations)
225 :
226 3497 : CALL strip_control_codes(csym%international_symbol)
227 3497 : CALL strip_control_codes(csym%schoenflies)
228 3497 : CALL strip_control_codes(csym%pointgroup_symbol)
229 : END IF
230 : END IF
231 : END IF
232 3808 : csym%symlib = spglib
233 :
234 3808 : CALL timestop(handle)
235 :
236 3808 : 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 3492 : 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 3492 : 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 3492 : CALL timeset(routineN, handle)
271 :
272 3492 : IF (PRESENT(shift)) THEN
273 13968 : csym%wvk0 = shift
274 : ELSE
275 0 : csym%wvk0 = 0.0_dp
276 : END IF
277 :
278 3492 : csym%istriz = -1
279 3492 : IF (PRESENT(symm)) THEN
280 3492 : IF (symm) csym%istriz = 1
281 : END IF
282 :
283 3492 : IF (PRESENT(full_grid)) THEN
284 3480 : fullmesh = full_grid
285 : ELSE
286 : fullmesh = .FALSE.
287 : END IF
288 3492 : csym%fullgrid = fullmesh
289 :
290 3492 : IF (PRESENT(gamma_centered)) THEN
291 3480 : gamma_mesh = gamma_centered
292 : ELSE
293 12 : gamma_mesh = .FALSE.
294 : END IF
295 :
296 3492 : IF (PRESENT(inversion_symmetry_only)) THEN
297 3480 : inversion_only = inversion_symmetry_only
298 : ELSE
299 : inversion_only = .FALSE.
300 : END IF
301 3492 : csym%inversion_only = inversion_only
302 :
303 3492 : IF (PRESENT(use_spglib_reduction)) THEN
304 3480 : spglib_reduction = use_spglib_reduction
305 : ELSE
306 12 : spglib_reduction = .FALSE.
307 : END IF
308 3492 : csym%spglib_reduction = spglib_reduction
309 :
310 3492 : IF (PRESENT(use_spglib_backend)) THEN
311 3480 : spglib_backend = use_spglib_backend
312 : ELSE
313 : spglib_backend = .FALSE.
314 : END IF
315 3492 : csym%spglib_backend = spglib_backend
316 :
317 3492 : 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 3492 : (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 3492 : csym%nkpoint = 0
328 13968 : csym%mesh(1:3) = nk(1:3)
329 3492 : csym%nrtot = 0
330 3492 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
331 3492 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
332 3492 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
333 3492 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
334 6984 : ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
335 :
336 3492 : nkpts = nk(1)*nk(2)*nk(3)
337 24444 : ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
338 : ! kp: link
339 10476 : ALLOCATE (csym%kplink(2, nkpts))
340 108120 : csym%kplink = 0
341 3492 : kpop = 0
342 :
343 : ! go through all the options
344 3492 : IF (csym%symlib) THEN
345 : ! symmetry library is available
346 3186 : IF (fullmesh) THEN
347 : ! full mesh requested
348 348 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
349 348 : IF (csym%istriz == 1) THEN
350 : ! use inversion symmetry
351 348 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :), nk)
352 : ELSE
353 : ! full kpoint mesh is used
354 : END IF
355 2838 : ELSE IF (csym%istriz /= 1 .OR. inversion_only) THEN
356 : ! use inversion symmetry
357 1104 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
358 1104 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :), nk)
359 : ELSE
360 : ! use symmetry library to reduce k-points
361 1734 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
362 1734 : IF (spglib_backend) THEN
363 760 : CALL kp_symmetry_spglib(csym, xkp, wkp, kpop)
364 : ELSE
365 974 : 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 306 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
372 306 : 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 294 : ELSE IF (csym%istriz /= 1 .AND. fullmesh) THEN
376 : ! full kpoint mesh is used
377 1242 : DO i = 1, nkpts
378 1242 : csym%kplink(1, i) = i
379 : END DO
380 : ELSE
381 : ! use inversion symmetry
382 164 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :), nk)
383 : END IF
384 : END IF
385 : ! count kpoints
386 3492 : nkp = 0
387 38368 : DO i = 1, nkpts
388 38368 : IF (wkp(i) > 0.0_dp) nkp = nkp + 1
389 : END DO
390 :
391 : ! store reduced kpoint set
392 3492 : csym%nkpoint = nkp
393 17460 : ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
394 6984 : ALLOCATE (rep_to_red(nkpts))
395 3492 : rep_to_red(:) = 0
396 3492 : j = 0
397 38368 : DO ik = 1, nkpts
398 38368 : IF (wkp(ik) > 0.0_dp) THEN
399 12462 : j = j + 1
400 12462 : csym%wkpoint(j) = wkp(ik)
401 49848 : csym%xkpoint(1:3, j) = xkp(1:3, ik)
402 12462 : rep_to_red(ik) = j
403 : END IF
404 : END DO
405 3492 : CPASSERT(j == nkp)
406 :
407 : ! kp: mesh
408 6984 : ALLOCATE (csym%kpmesh(3, nkpts))
409 142996 : csym%kpmesh(1:3, 1:nkpts) = xkp(1:3, 1:nkpts)
410 :
411 : ! kp: link
412 38368 : DO ik = 1, nkpts
413 34876 : i = csym%kplink(1, ik)
414 34876 : CPASSERT(i > 0 .AND. i <= nkpts)
415 34876 : csym%kplink(2, ik) = rep_to_red(i)
416 38368 : CPASSERT(csym%kplink(2, ik) > 0)
417 : END DO
418 3492 : DEALLOCATE (rep_to_red)
419 :
420 : ! kp: operations
421 6984 : ALLOCATE (csym%kpop(nkpts))
422 3492 : IF (csym%nrtot > 0 .AND. csym%istriz == 1 .AND. .NOT. fullmesh .AND. &
423 : .NOT. inversion_only) THEN
424 : ! atomic symmetry operations possible
425 21900 : csym%kpop(1:nkpts) = kpop(1:nkpts)
426 21900 : DO ik = 1, nkpts
427 21900 : CPASSERT(csym%kpop(ik) /= 0)
428 : END DO
429 : ELSE
430 : ! only time reversal symmetry
431 16468 : DO ik = 1, nkpts
432 16468 : IF (wkp(ik) > 0.0_dp) THEN
433 8204 : csym%kpop(ik) = 1
434 : ELSE
435 6516 : csym%kpop(ik) = 2
436 : END IF
437 : END DO
438 : END IF
439 :
440 3492 : DEALLOCATE (xkp, wkp, kpop)
441 :
442 3492 : CALL timestop(handle)
443 :
444 3492 : 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 986 : 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 986 : 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 986 : IF (PRESENT(use_spglib_reduction)) THEN
580 986 : spglib_reduction = use_spglib_reduction
581 : ELSE
582 : spglib_reduction = .FALSE.
583 : END IF
584 :
585 986 : CALL detect_k290_symmetry(csym, a1, a2, a3, b1, b2, b3, alat, nc, ib, r)
586 986 : IF (spglib_reduction) CALL cite_reference(Togo2018)
587 :
588 986 : IF (csym%n_operations > nc .AND. .NOT. spglib_reduction) THEN
589 714 : ALLOCATE (srot(3, 3, csym%n_operations))
590 238 : CALL setup_spglib_operations(csym, srot, nrot)
591 238 : CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
592 476 : DEALLOCATE (srot)
593 740 : 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 732 : CALL reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
601 : END IF
602 :
603 986 : 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 994 : 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 994 : INTEGER, ALLOCATABLE, DIMENSION(:) :: isc, ty
686 994 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0
687 994 : 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 994 : iou = csym%punit
692 994 : nat = csym%nat
693 994 : CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
694 :
695 9940 : ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, nat), ty(nat), isc(nat), f0(49, nat))
696 6176 : ty(1:nat) = csym%atype(1:nat)
697 6176 : DO i = 1, nat
698 89088 : xkapa(1:3, i) = MATMUL(csym%hmat, csym%scoord(1:3, i))
699 : END DO
700 :
701 994 : IF (iou > 0) THEN
702 : WRITE (iou, '(/,(T2,A79))') &
703 340 : "*******************************************************************************", &
704 340 : "** Crystal Symmetry Analysis by K290 **", &
705 680 : "*******************************************************************************"
706 : END IF
707 994 : 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 994 : vt, f0, r, tvec, origin, rx, isc, csym%delta)
712 :
713 994 : IF (iou > 0) THEN
714 : WRITE (iou, '((T2,A79))') &
715 340 : "*******************************************************************************", &
716 340 : "** Finished K290 **", &
717 680 : "*******************************************************************************"
718 : END IF
719 :
720 994 : csym%nrtot = nc
721 994 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
722 994 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
723 994 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
724 994 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
725 9940 : ALLOCATE (csym%rt(3, 3, nc), csym%vt(3, nc), csym%ibrot(nc), csym%f0(nat, nc))
726 68354 : csym%vt(1:3, 1:nc) = vt(1:3, 1:nc)
727 17834 : DO i = 1, nc
728 218920 : csym%rt(1:3, 1:3, i) = r(1:3, 1:3, ib(i))
729 93618 : csym%f0(1:nat, i) = f0(i, 1:nat)
730 : END DO
731 17834 : csym%ibrot(1:nc) = ib(1:nc)
732 :
733 994 : DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
734 :
735 994 : 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 1002 : 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 4008 : a1(1:3) = csym%hmat(1:3, 1)
756 4008 : a2(1:3) = csym%hmat(1:3, 2)
757 4008 : a3(1:3) = csym%hmat(1:3, 3)
758 4008 : 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 1002 : a2(3)*a3(2)*a1(1) - a3(3)*a1(2)*a2(1)
762 1002 : volum = ABS(volum)
763 1002 : b1(1) = (a2(2)*a3(3) - a2(3)*a3(2))/volum
764 1002 : b1(2) = (a2(3)*a3(1) - a2(1)*a3(3))/volum
765 1002 : b1(3) = (a2(1)*a3(2) - a2(2)*a3(1))/volum
766 1002 : b2(1) = (a3(2)*a1(3) - a3(3)*a1(2))/volum
767 1002 : b2(2) = (a3(3)*a1(1) - a3(1)*a1(3))/volum
768 1002 : b2(3) = (a3(1)*a1(2) - a3(2)*a1(1))/volum
769 1002 : b3(1) = (a1(2)*a2(3) - a1(3)*a2(2))/volum
770 1002 : b3(2) = (a1(3)*a2(1) - a1(1)*a2(3))/volum
771 1002 : b3(3) = (a1(1)*a2(2) - a1(2)*a2(1))/volum
772 :
773 1002 : 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 1016 : 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 1016 : 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 1016 : CPASSERT(csym%symlib)
795 :
796 653564 : srot = 0
797 1016 : csym%nrtot = 0
798 1016 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
799 1016 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
800 1016 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
801 1016 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
802 5080 : ALLOCATE (csym%rt(3, 3, csym%n_operations), csym%vt(3, csym%n_operations))
803 6096 : ALLOCATE (csym%ibrot(csym%n_operations), csym%f0(csym%nat, csym%n_operations))
804 653564 : csym%rt = 0.0_dp
805 201800 : csym%vt = 0.0_dp
806 51212 : csym%ibrot = 0
807 417836 : csym%f0 = 0
808 :
809 1016 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
810 1016 : h_inv = inv_3x3(csym%hmat)
811 3048 : ALLOCATE (perm(csym%nat))
812 :
813 1016 : eye = 0
814 1016 : eye(1, 1) = 1
815 1016 : eye(2, 2) = 1
816 1016 : eye(3, 3) = 1
817 :
818 1016 : 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 5080 : DO pass = 1, 4
822 205864 : DO iop = 1, csym%n_operations
823 2610192 : irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
824 2610192 : frot(1:3, 1:3) = TRANSPOSE(irot(1:3, 1:3))
825 458992 : identity = ALL(frot == eye)
826 : zero_translation = ALL(ABS(csym%translations(1:3, iop) - &
827 348184 : ANINT(csym%translations(1:3, iop))) <= eps)
828 200784 : IF (pass == 1 .AND. (.NOT. identity .OR. .NOT. zero_translation)) CYCLE
829 151604 : IF (pass == 2 .AND. (identity .OR. .NOT. zero_translation)) CYCLE
830 110814 : IF (pass == 3 .AND. (.NOT. identity .OR. zero_translation)) CYCLE
831 61418 : IF (pass == 4 .AND. (identity .OR. zero_translation)) CYCLE
832 :
833 50196 : duplicate = .FALSE.
834 1143540 : DO jop = 1, nrot
835 2790048 : IF (ALL(frot == srot(:, :, jop))) THEN
836 : duplicate = .TRUE.
837 : EXIT
838 : END IF
839 : END DO
840 50196 : IF (duplicate) CYCLE
841 :
842 16884 : CALL spglib_atom_permutation(csym, frot, csym%translations(:, iop), perm, valid)
843 16884 : IF (.NOT. valid) CYCLE
844 :
845 16884 : nrot = nrot + 1
846 :
847 219492 : srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
848 219492 : rfrac(1:3, 1:3) = REAL(frot(1:3, 1:3), KIND=dp)
849 1536444 : csym%rt(1:3, 1:3, nrot) = MATMUL(csym%hmat, MATMUL(rfrac, h_inv))
850 67536 : csym%vt(1:3, nrot) = csym%translations(1:3, iop)
851 16884 : csym%ibrot(nrot) = nrot
852 123380 : csym%f0(1:csym%nat, nrot) = perm(1:csym%nat)
853 : END DO
854 : END DO
855 :
856 1016 : DEALLOCATE (perm)
857 1016 : csym%nrtot = nrot
858 1016 : IF (nrot == 0) CALL cp_abort(__LOCATION__, "SPGLIB did not return usable symmetry operations")
859 :
860 1016 : 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 16884 : 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 16884 : 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 16884 : nat = csym%nat
936 16884 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
937 219492 : rfrac(1:3, 1:3) = REAL(rot(1:3, 1:3), KIND=dp)
938 50652 : ALLOCATE (used(nat))
939 16884 : used = .FALSE.
940 119316 : perm = 0
941 16884 : valid = .TRUE.
942 :
943 119316 : DO i = 1, nat
944 1638912 : spos(1:3) = MATMUL(rfrac(1:3, 1:3), csym%scoord(1:3, i)) + trans(1:3)
945 459424 : found = .FALSE.
946 459424 : DO j = 1, nat
947 459424 : IF (used(j)) CYCLE
948 251330 : IF (csym%atype(i) /= csym%atype(j)) CYCLE
949 1005320 : diff(1:3) = spos(1:3) - csym%scoord(1:3, j)
950 1005320 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
951 577650 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
952 102432 : perm(i) = j
953 102432 : used(j) = .TRUE.
954 : found = .TRUE.
955 : EXIT
956 : END IF
957 : END DO
958 16884 : IF (.NOT. found) THEN
959 0 : valid = .FALSE.
960 0 : EXIT
961 : END IF
962 : END DO
963 :
964 16884 : DEALLOCATE (used)
965 :
966 16884 : 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 998 : 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 998 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kscore, op_score
987 998 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: krot
988 : REAL(KIND=dp) :: eps
989 : REAL(KIND=dp), DIMENSION(3) :: rr
990 :
991 998 : nkpts = SIZE(wkp)
992 998 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
993 6986 : ALLOCATE (kscore(nkpts), op_score(nrot), krot(3, 3, nrot))
994 :
995 17018 : DO iop = 1, nrot
996 16020 : kr = csym%ibrot(iop)
997 208260 : krot(1:3, 1:3, iop) = reciprocal_rotation(srot(1:3, 1:3, kr))
998 17018 : op_score(iop) = spglib_operation_score(csym, iop, srot(1:3, 1:3, kr))
999 : END DO
1000 :
1001 14072 : wkp = 0.0_dp
1002 14072 : kpop = 0
1003 14072 : csym%kplink(1, :) = 0
1004 14072 : kscore = HUGE(0)
1005 :
1006 14072 : DO i = 1, nkpts
1007 13074 : IF (csym%kplink(1, i) /= 0) CYCLE
1008 :
1009 2442 : csym%kplink(1, i) = i
1010 2442 : wkp(i) = 1.0_dp
1011 2442 : kpop(i) = 1
1012 2442 : kscore(i) = 0
1013 :
1014 43714 : DO iop = 1, nrot
1015 133896 : DO isign = 1, 2
1016 2013700 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3, iop), KIND=dp), xkp(1:3, i))
1017 80548 : kr = csym%ibrot(iop)
1018 80548 : IF (isign == 2) THEN
1019 161096 : rr(1:3) = -rr(1:3)
1020 40274 : kr = -kr
1021 : END IF
1022 :
1023 80548 : j = regular_kpoint_index(rr, xkp(1:3, 1), csym%mesh, eps)
1024 80548 : IF (j == 0) CYCLE
1025 113782 : IF (csym%kplink(1, j) == 0) THEN
1026 10632 : csym%kplink(1, j) = i
1027 10632 : wkp(i) = wkp(i) + 1.0_dp
1028 10632 : kpop(j) = kr
1029 10632 : kscore(j) = op_score(iop)
1030 : ELSE
1031 62876 : CPASSERT(csym%kplink(1, j) == i)
1032 62876 : IF (op_score(iop) < kscore(j)) THEN
1033 3524 : kpop(j) = kr
1034 3524 : kscore(j) = op_score(iop)
1035 : END IF
1036 : END IF
1037 : END DO
1038 : END DO
1039 : END DO
1040 :
1041 14072 : DO i = 1, nkpts
1042 13074 : CPASSERT(csym%kplink(1, i) /= 0)
1043 14072 : CPASSERT(kpop(i) /= 0)
1044 : END DO
1045 998 : DEALLOCATE (kscore, op_score, krot)
1046 :
1047 998 : 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 16020 : 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 16020 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1580 16020 : nat = SIZE(csym%f0, 1)
1581 16020 : score = 0
1582 111540 : DO i = 1, nat
1583 111540 : IF (csym%f0(i, iop) /= i) score = score + 100
1584 : END DO
1585 46114 : IF (ANY(ABS(csym%vt(1:3, iop) - ANINT(csym%vt(1:3, iop))) > eps)) score = score + 10
1586 :
1587 16020 : eye = 0
1588 16020 : eye(1, 1) = 1
1589 16020 : eye(2, 2) = 1
1590 16020 : eye(3, 3) = 1
1591 640800 : r2(1:3, 1:3) = MATMUL(srot(1:3, 1:3), srot(1:3, 1:3))
1592 113416 : IF (ANY(r2(1:3, 1:3) /= eye(1:3, 1:3))) score = score + 1
1593 :
1594 16020 : 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 17684 : 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 229892 : rinv = inv_3x3(REAL(rot(1:3, 1:3), KIND=dp))
1608 229892 : krot(1:3, 1:3) = NINT(TRANSPOSE(rinv(1:3, 1:3)))
1609 :
1610 17684 : 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 732 : 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 732 : nkpts = SIZE(wkp)
1646 732 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1647 732 : nskipped_overlaps = 0
1648 :
1649 7516 : wkp = 0.0_dp
1650 7516 : kpop = 0
1651 7516 : csym%kplink(1, :) = 0
1652 :
1653 7516 : DO i = 1, nkpts
1654 6784 : IF (csym%kplink(1, i) /= 0) CYCLE
1655 :
1656 1766 : csym%kplink(1, i) = i
1657 1766 : wkp(i) = 1.0_dp
1658 1766 : kpop(i) = 1
1659 :
1660 27510 : DO iop = 1, nc
1661 81820 : DO ibsign = 1, 2
1662 50024 : kr = ib(iop)
1663 200096 : wcart(1:3) = alat*(xkp(1, i)*b1(1:3) + xkp(2, i)*b2(1:3) + xkp(3, i)*b3(1:3))
1664 200096 : wcart(1:3) = kp_apply_operation(wcart(1:3), r(1:3, 1:3, kr))
1665 50024 : IF (ibsign == 2) THEN
1666 100048 : wcart(1:3) = -wcart(1:3)
1667 25012 : kr = -kr
1668 : END IF
1669 200096 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1670 200096 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1671 200096 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1672 :
1673 50024 : 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 50024 : IF (j == 0) CYCLE
1676 74192 : IF (csym%kplink(1, j) == 0) THEN
1677 5018 : csym%kplink(1, j) = i
1678 5018 : wkp(i) = wkp(i) + 1.0_dp
1679 5018 : kpop(j) = kr
1680 44162 : 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 732 : IF (nskipped_overlaps > 0) THEN
1691 56 : wkp = 1.0_dp
1692 2 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :), csym%mesh)
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 7460 : DO i = 1, nkpts
1708 6730 : CPASSERT(csym%kplink(1, i) /= 0)
1709 7460 : CPASSERT(kpop(i) /= 0)
1710 : END DO
1711 : END IF
1712 :
1713 732 : 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 145396 : 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 145396 : index = 0
1734 581584 : IF (ANY(mesh <= 0)) RETURN
1735 564132 : DO idim = 1, 3
1736 426620 : scaled = (x(idim) - origin(idim))*REAL(mesh(idim), KIND=dp)
1737 426620 : istep(idim) = NINT(scaled)
1738 426620 : IF (ABS(scaled - REAL(istep(idim), KIND=dp)) > eps*REAL(mesh(idim), KIND=dp)) RETURN
1739 556248 : icoord(idim) = MODULO(istep(idim), mesh(idim)) + 1
1740 : END DO
1741 :
1742 137512 : index = ((icoord(1) - 1)*mesh(2) + icoord(2) - 1)*mesh(3) + icoord(3)
1743 :
1744 137512 : 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 3492 : 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 3492 : IF (PRESENT(gamma_centered)) THEN
1767 3492 : gamma_mesh = gamma_centered
1768 : ELSE
1769 : gamma_mesh = .FALSE.
1770 : END IF
1771 :
1772 38368 : wkp = 0.0_dp
1773 3492 : i = 0
1774 10406 : DO ix = 1, nk(1)
1775 24978 : DO iy = 1, nk(2)
1776 56362 : DO iz = 1, nk(3)
1777 34876 : i = i + 1
1778 34876 : ik(1) = ix
1779 34876 : ik(2) = iy
1780 34876 : ik(3) = iz
1781 139504 : DO idim = 1, 3
1782 139504 : 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 102480 : (2._dp*REAL(nk(idim), KIND=dp))
1788 : END IF
1789 : END DO
1790 139504 : xkp(1:3, i) = kpt_latt(1:3)
1791 49448 : wkp(i) = 1.0_dp
1792 : END DO
1793 : END DO
1794 : END DO
1795 38368 : DO i = 1, nk(1)*nk(2)*nk(3)
1796 142996 : xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
1797 : END DO
1798 :
1799 3492 : END SUBROUTINE full_grid_gen
1800 :
1801 : ! **************************************************************************************************
1802 : !> \brief Reduce a regular mesh by inversion/time-reversal in linear time.
1803 : !> \param xkp full mesh in the order produced by full_grid_gen
1804 : !> \param wkp reduced k-point weights
1805 : !> \param link full-mesh index of each point's representative
1806 : !> \param mesh mesh dimensions
1807 : ! **************************************************************************************************
1808 1618 : SUBROUTINE inversion_symm(xkp, wkp, link, mesh)
1809 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1810 : REAL(KIND=dp), DIMENSION(:) :: wkp
1811 : INTEGER, DIMENSION(:) :: link
1812 : INTEGER, DIMENSION(3), INTENT(IN) :: mesh
1813 :
1814 : INTEGER :: i, j, nkpts
1815 : REAL(KIND=dp), DIMENSION(3) :: diff
1816 :
1817 1618 : nkpts = SIZE(wkp, 1)
1818 :
1819 15226 : link(:) = 0
1820 15226 : DO i = 1, nkpts
1821 13608 : IF (link(i) == 0) link(i) = i
1822 54432 : j = regular_kpoint_index(-xkp(1:3, i), xkp(1:3, 1), mesh, 2.e-12_dp)
1823 : ! Keep the first point as representative; zero means -k is off the shifted mesh.
1824 13608 : IF (j <= i) CYCLE
1825 6516 : IF (wkp(j) == 0.0_dp) CYCLE
1826 : ! Use the original coordinate tolerance after allowing for rounding in the index lookup.
1827 26064 : diff(1:3) = xkp(1:3, i) + xkp(1:3, j)
1828 26064 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1829 27682 : IF (ALL(ABS(diff(1:3)) < 1.e-12_dp)) THEN
1830 6516 : wkp(i) = wkp(i) + wkp(j)
1831 6516 : wkp(j) = 0.0_dp
1832 6516 : link(j) = i
1833 : END IF
1834 : END DO
1835 :
1836 1618 : END SUBROUTINE inversion_symm
1837 :
1838 : ! **************************************************************************************************
1839 : !> \brief ...
1840 : !> \param x ...
1841 : !> \param r ...
1842 : !> \return ...
1843 : ! **************************************************************************************************
1844 53336 : FUNCTION kp_apply_operation(x, r) RESULT(y)
1845 : REAL(KIND=dp), INTENT(IN) :: x(3), r(3, 3)
1846 : REAL(KIND=dp) :: y(3)
1847 :
1848 53336 : y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
1849 53336 : y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
1850 53336 : y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
1851 :
1852 53336 : END FUNCTION kp_apply_operation
1853 :
1854 : ! **************************************************************************************************
1855 : !> \brief ...
1856 : !> \param csym ...
1857 : ! **************************************************************************************************
1858 3502 : SUBROUTINE print_crys_symmetry(csym)
1859 : TYPE(csym_type) :: csym
1860 :
1861 : INTEGER :: i, iunit, j, plevel
1862 :
1863 3502 : iunit = csym%punit
1864 3502 : IF (iunit >= 0) THEN
1865 1473 : plevel = csym%plevel
1866 1473 : WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
1867 1473 : IF (csym%symlib) THEN
1868 1471 : WRITE (iunit, '(A,T71,A10)') " International Symbol: ", ADJUSTR(TRIM(csym%international_symbol))
1869 1471 : WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", ADJUSTR(TRIM(csym%pointgroup_symbol))
1870 1471 : WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", ADJUSTR(TRIM(csym%schoenflies))
1871 : !
1872 1471 : WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
1873 1471 : IF (plevel > 0) THEN
1874 0 : DO i = 1, csym%n_operations
1875 : WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
1876 0 : " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
1877 0 : WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
1878 : END DO
1879 : END IF
1880 : ELSE
1881 2 : IF (csym%spglib_requested) THEN
1882 1 : WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not available"
1883 : ELSE
1884 1 : WRITE (iunit, "(T2,A)") "SPGLIB Crystal Symmetry Information was not requested"
1885 : END IF
1886 : END IF
1887 : END IF
1888 :
1889 3502 : END SUBROUTINE print_crys_symmetry
1890 :
1891 : ! **************************************************************************************************
1892 : !> \brief ...
1893 : !> \param csym ...
1894 : ! **************************************************************************************************
1895 3212 : SUBROUTINE print_kp_symmetry(csym)
1896 : TYPE(csym_type), INTENT(IN) :: csym
1897 :
1898 : INTEGER :: i, iunit, nat, nmesh, plevel
1899 :
1900 3212 : iunit = csym%punit
1901 3212 : IF (iunit >= 0) THEN
1902 1183 : plevel = csym%plevel
1903 1183 : WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
1904 1183 : WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
1905 1183 : WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
1906 4584 : DO i = 1, csym%nkpoint
1907 4584 : WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), NINT(csym%wkpoint(i))
1908 : END DO
1909 1183 : nmesh = csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
1910 1183 : IF (nmesh > 0) THEN
1911 1175 : WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
1912 : ELSE
1913 8 : nmesh = SIZE(csym%kpmesh, 2)
1914 8 : WRITE (iunit, '(/,A,T70,I10)') " Explicit K-point Set: ", nmesh
1915 : END IF
1916 1183 : WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
1917 11331 : DO i = 1, nmesh
1918 10148 : WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
1919 21479 : csym%kplink(1:2, i), csym%kpop(i)
1920 : END DO
1921 1183 : IF (csym%nrtot > 0) THEN
1922 667 : WRITE (iunit, '(/,A)') " Atom Transformation Table"
1923 667 : nat = SIZE(csym%f0, 1)
1924 7692 : DO i = 1, csym%nrtot
1925 7692 : WRITE (iunit, '(T10,A,I5,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
1926 : END DO
1927 : END IF
1928 : END IF
1929 :
1930 3212 : END SUBROUTINE print_kp_symmetry
1931 :
1932 0 : END MODULE cryssym
|
