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