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 3491 : SUBROUTINE release_csym_type(csym)
90 : TYPE(csym_type) :: csym
91 :
92 3491 : IF (ALLOCATED(csym%rotations)) THEN
93 3314 : DEALLOCATE (csym%rotations)
94 : END IF
95 3491 : IF (ALLOCATED(csym%translations)) THEN
96 3314 : DEALLOCATE (csym%translations)
97 : END IF
98 3491 : IF (ALLOCATED(csym%atype)) THEN
99 3491 : DEALLOCATE (csym%atype)
100 : END IF
101 3491 : IF (ALLOCATED(csym%scoord)) THEN
102 3491 : DEALLOCATE (csym%scoord)
103 : END IF
104 3491 : IF (ALLOCATED(csym%xkpoint)) THEN
105 3206 : DEALLOCATE (csym%xkpoint)
106 : END IF
107 3491 : IF (ALLOCATED(csym%wkpoint)) THEN
108 3206 : DEALLOCATE (csym%wkpoint)
109 : END IF
110 3491 : IF (ALLOCATED(csym%kpmesh)) THEN
111 3206 : DEALLOCATE (csym%kpmesh)
112 : END IF
113 3491 : IF (ALLOCATED(csym%kplink)) THEN
114 3206 : DEALLOCATE (csym%kplink)
115 : END IF
116 3491 : IF (ALLOCATED(csym%kpop)) THEN
117 3206 : DEALLOCATE (csym%kpop)
118 : END IF
119 3491 : IF (ALLOCATED(csym%rt)) THEN
120 3206 : DEALLOCATE (csym%rt)
121 : END IF
122 3491 : IF (ALLOCATED(csym%vt)) THEN
123 3206 : DEALLOCATE (csym%vt)
124 : END IF
125 3491 : IF (ALLOCATED(csym%f0)) THEN
126 3206 : DEALLOCATE (csym%f0)
127 : END IF
128 3491 : IF (ALLOCATED(csym%ibrot)) THEN
129 3206 : DEALLOCATE (csym%ibrot)
130 : END IF
131 :
132 3491 : 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 3491 : 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 3491 : CALL timeset(routineN, handle)
160 :
161 : !..total number of atoms
162 3491 : nat = SIZE(scoor, 2)
163 3491 : csym%nat = nat
164 :
165 : ! output unit
166 3491 : IF (PRESENT(iounit)) THEN
167 3491 : csym%punit = iounit
168 : ELSE
169 0 : csym%punit = -1
170 : END IF
171 :
172 : ! accuracy for symmetry
173 3491 : IF (PRESENT(delta)) THEN
174 3491 : csym%delta = delta
175 : ELSE
176 0 : csym%delta = 1.e-6_dp
177 : END IF
178 :
179 : !..set cell values
180 45383 : csym%hmat = hmat
181 :
182 : ! atom types
183 10473 : ALLOCATE (csym%atype(nat))
184 23100 : csym%atype(1:nat) = types(1:nat)
185 :
186 : ! scaled coordinates
187 10473 : ALLOCATE (csym%scoord(3, nat))
188 81927 : csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
189 :
190 3491 : csym%n_operations = 0
191 :
192 : !..try spglib
193 3491 : my_use_spglib = .TRUE.
194 3491 : IF (PRESENT(use_spglib)) my_use_spglib = use_spglib
195 285 : csym%spglib_requested = my_use_spglib
196 3206 : IF (.NOT. my_use_spglib) THEN
197 : spglib = .FALSE.
198 : ELSE
199 3315 : major = spg_get_major_version()
200 3315 : minor = spg_get_minor_version()
201 3315 : micro = spg_get_micro_version()
202 3315 : IF (major == 0) THEN
203 0 : CALL cp_warn(__LOCATION__, "Symmetry library SPGLIB not available")
204 0 : spglib = .FALSE.
205 : ELSE
206 3315 : spglib = .TRUE.
207 3315 : CALL cite_reference(Togo2018)
208 3315 : ierr = spg_get_international(csym%international_symbol, TRANSPOSE(hmat), scoor, types, nat, delta)
209 3315 : IF (ierr == 0) THEN
210 1 : CALL cp_warn(__LOCATION__, "Symmetry Library SPGLIB failed")
211 1 : spglib = .FALSE.
212 : ELSE
213 3314 : nop = spg_get_multiplicity(TRANSPOSE(hmat), scoor, types, nat, delta)
214 16570 : ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
215 3314 : csym%n_operations = nop
216 : ierr = spg_get_symmetry(csym%rotations, csym%translations, nop, &
217 3314 : TRANSPOSE(hmat), scoor, types, nat, delta)
218 : ! Schoenflies Symbol
219 3314 : csym%schoenflies = ' '
220 3314 : ierr = spg_get_schoenflies(csym%schoenflies, TRANSPOSE(hmat), scoor, types, nat, delta)
221 : ! Point Group
222 3314 : csym%pointgroup_symbol = ' '
223 3314 : tra_mat = 0
224 : ierr = spg_get_pointgroup(csym%pointgroup_symbol, tra_mat, &
225 3314 : csym%rotations, csym%n_operations)
226 :
227 3314 : CALL strip_control_codes(csym%international_symbol)
228 3314 : CALL strip_control_codes(csym%schoenflies)
229 3314 : CALL strip_control_codes(csym%pointgroup_symbol)
230 : END IF
231 : END IF
232 : END IF
233 3491 : csym%symlib = spglib
234 :
235 3491 : CALL timestop(handle)
236 :
237 3491 : 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 3180 : 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 3180 : 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 3180 : CALL timeset(routineN, handle)
272 :
273 3180 : IF (PRESENT(shift)) THEN
274 12720 : csym%wvk0 = shift
275 : ELSE
276 0 : csym%wvk0 = 0.0_dp
277 : END IF
278 :
279 3180 : csym%istriz = -1
280 3180 : IF (PRESENT(symm)) THEN
281 3180 : IF (symm) csym%istriz = 1
282 : END IF
283 :
284 3180 : IF (PRESENT(full_grid)) THEN
285 3180 : fullmesh = full_grid
286 : ELSE
287 : fullmesh = .FALSE.
288 : END IF
289 3180 : csym%fullgrid = fullmesh
290 :
291 3180 : IF (PRESENT(gamma_centered)) THEN
292 3180 : gamma_mesh = gamma_centered
293 : ELSE
294 0 : gamma_mesh = .FALSE.
295 : END IF
296 :
297 3180 : IF (PRESENT(inversion_symmetry_only)) THEN
298 3180 : inversion_only = inversion_symmetry_only
299 : ELSE
300 : inversion_only = .FALSE.
301 : END IF
302 3180 : csym%inversion_only = inversion_only
303 :
304 3180 : IF (PRESENT(use_spglib_reduction)) THEN
305 3180 : spglib_reduction = use_spglib_reduction
306 : ELSE
307 0 : spglib_reduction = .FALSE.
308 : END IF
309 3180 : csym%spglib_reduction = spglib_reduction
310 :
311 3180 : IF (PRESENT(use_spglib_backend)) THEN
312 3180 : spglib_backend = use_spglib_backend
313 : ELSE
314 : spglib_backend = .FALSE.
315 : END IF
316 3180 : csym%spglib_backend = spglib_backend
317 :
318 3180 : 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 3180 : (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 3180 : csym%nkpoint = 0
329 12720 : csym%mesh(1:3) = nk(1:3)
330 3180 : csym%nrtot = 0
331 3180 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
332 3180 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
333 3180 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
334 3180 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
335 6360 : ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
336 :
337 3180 : nkpts = nk(1)*nk(2)*nk(3)
338 22260 : ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
339 : ! kp: link
340 9540 : ALLOCATE (csym%kplink(2, nkpts))
341 96234 : csym%kplink = 0
342 3180 : kpop = 0
343 :
344 : ! go through all the options
345 3180 : IF (csym%symlib) THEN
346 : ! symmetry library is available
347 3008 : IF (fullmesh) THEN
348 : ! full mesh requested
349 342 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
350 342 : IF (csym%istriz == 1) THEN
351 : ! use inversion symmetry
352 342 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
353 : ELSE
354 : ! full kpoint mesh is used
355 : END IF
356 2666 : ELSE IF (csym%istriz /= 1 .OR. inversion_only) THEN
357 : ! use inversion symmetry
358 1226 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
359 1226 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
360 : ELSE
361 : ! use symmetry library to reduce k-points
362 1440 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
363 1440 : IF (spglib_backend) THEN
364 720 : CALL kp_symmetry_spglib(csym, xkp, wkp, kpop)
365 : ELSE
366 720 : 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 172 : CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
373 172 : 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 172 : ELSE IF (csym%istriz /= 1 .AND. fullmesh) THEN
377 : ! full kpoint mesh is used
378 870 : DO i = 1, nkpts
379 870 : csym%kplink(1, i) = i
380 : END DO
381 : ELSE
382 : ! use inversion symmetry
383 102 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
384 : END IF
385 : END IF
386 : ! count kpoints
387 3180 : nkp = 0
388 34198 : DO i = 1, nkpts
389 34198 : IF (wkp(i) > 0.0_dp) nkp = nkp + 1
390 : END DO
391 :
392 : ! store reduced kpoint set
393 3180 : csym%nkpoint = nkp
394 15900 : ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
395 9540 : ALLOCATE (xptr(nkp))
396 34198 : j = 0
397 34198 : DO ik = 1, nkpts
398 34198 : IF (wkp(ik) > 0.0_dp) THEN
399 11442 : j = j + 1
400 11442 : csym%wkpoint(j) = wkp(ik)
401 45768 : csym%xkpoint(1:3, j) = xkp(1:3, ik)
402 11442 : xptr(j) = ik
403 : END IF
404 : END DO
405 3180 : CPASSERT(j == nkp)
406 :
407 : ! kp: mesh
408 6360 : ALLOCATE (csym%kpmesh(3, nkpts))
409 127252 : csym%kpmesh(1:3, 1:nkpts) = xkp(1:3, 1:nkpts)
410 :
411 : ! kp: link
412 34198 : DO ik = 1, nkpts
413 31018 : i = csym%kplink(1, ik)
414 177630 : DO j = 1, nkp
415 174450 : IF (i == xptr(j)) THEN
416 31018 : csym%kplink(2, ik) = j
417 31018 : EXIT
418 : END IF
419 : END DO
420 : END DO
421 3180 : DEALLOCATE (xptr)
422 :
423 : ! kp: operations
424 6360 : ALLOCATE (csym%kpop(nkpts))
425 3180 : IF (csym%nrtot > 0 .AND. csym%istriz == 1 .AND. .NOT. fullmesh .AND. &
426 : .NOT. inversion_only) THEN
427 : ! atomic symmetry operations possible
428 17470 : csym%kpop(1:nkpts) = kpop(1:nkpts)
429 17470 : DO ik = 1, nkpts
430 17470 : CPASSERT(csym%kpop(ik) /= 0)
431 : END DO
432 : ELSE
433 : ! only time reversal symmetry
434 16728 : DO ik = 1, nkpts
435 16728 : IF (wkp(ik) > 0.0_dp) THEN
436 8118 : csym%kpop(ik) = 1
437 : ELSE
438 6868 : csym%kpop(ik) = 2
439 : END IF
440 : END DO
441 : END IF
442 :
443 3180 : DEALLOCATE (xkp, wkp, kpop)
444 :
445 3180 : CALL timestop(handle)
446 :
447 3180 : 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 720 : 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 720 : INTEGER, ALLOCATABLE, DIMENSION(:) :: includ, isc, list, lwght, ty
578 720 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0, lrot
579 720 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
580 : INTEGER, DIMENSION(48) :: ib
581 : LOGICAL :: spglib_reduction
582 : REAL(KIND=dp) :: alat
583 720 : 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 720 : iou = csym%punit
590 9360 : hmat = csym%hmat
591 720 : nat = csym%nat
592 720 : iq1 = csym%mesh(1)
593 720 : iq2 = csym%mesh(2)
594 720 : iq3 = csym%mesh(3)
595 720 : 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 720 : wvk0 = 0.0_dp
599 720 : istriz = csym%istriz
600 720 : IF (PRESENT(use_spglib_reduction)) THEN
601 720 : spglib_reduction = use_spglib_reduction
602 : ELSE
603 : spglib_reduction = .FALSE.
604 : END IF
605 2880 : a1(1:3) = hmat(1:3, 1)
606 2880 : a2(1:3) = hmat(1:3, 2)
607 2880 : a3(1:3) = hmat(1:3, 3)
608 2880 : alat = SQRT(SUM(a1**2))
609 720 : strain = 0.0_dp
610 6480 : ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, 200), ty(nat), isc(nat), f0(49, nat))
611 4868 : ty(1:nat) = csym%atype(1:nat)
612 4868 : nsp = MAXVAL(ty)
613 4868 : DO i = 1, nat
614 67088 : xkapa(1:3, i) = MATMUL(hmat, csym%scoord(1:3, i))
615 : END DO
616 720 : nhash = MAX(1000, nkpoint)
617 5760 : ALLOCATE (wvkl(3, nkpoint), rlist(3, nkpoint), includ(nkpoint), list(nhash + nkpoint))
618 2880 : ALLOCATE (lrot(48, nkpoint), lwght(nkpoint))
619 :
620 720 : IF (iou > 0) THEN
621 : WRITE (iou, '(/,(T2,A79))') &
622 290 : "*******************************************************************************", &
623 290 : "** Special K-Point Generation by K290 **", &
624 580 : "*******************************************************************************"
625 : END IF
626 720 : CALL cite_reference(Worlton1972)
627 720 : 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 720 : 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 720 : vt, f0, r, tvec, origin, rx, isc, csym%delta)
637 :
638 720 : IF (iou > 0) THEN
639 : WRITE (iou, '((T2,A79))') &
640 290 : "*******************************************************************************", &
641 290 : "** Finished K290 **", &
642 580 : "*******************************************************************************"
643 : END IF
644 :
645 720 : csym%nrtot = nc
646 720 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
647 720 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
648 720 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
649 720 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
650 5040 : ALLOCATE (csym%rt(3, 3, nc), csym%vt(3, nc), csym%ibrot(nc))
651 34064 : csym%vt(1:3, 1:nc) = vt(1:3, 1:nc)
652 2880 : ALLOCATE (csym%f0(nat, nc))
653 9056 : DO i = 1, nc
654 108368 : csym%rt(1:3, 1:3, i) = r(1:3, 1:3, ib(i))
655 63880 : csym%f0(1:nat, i) = f0(i, 1:nat)
656 : END DO
657 9056 : csym%ibrot(1:nc) = ib(1:nc)
658 :
659 720 : IF (csym%n_operations > nc .AND. .NOT. spglib_reduction) THEN
660 : IF (ALLOCATED(srot)) DEALLOCATE (srot)
661 486 : ALLOCATE (srot(3, 3, csym%n_operations))
662 162 : CALL setup_spglib_operations(csym, srot, nrot)
663 162 : CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
664 324 : DEALLOCATE (srot)
665 550 : 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 542 : CALL reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
673 : END IF
674 720 : DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
675 720 : DEALLOCATE (wvkl, rlist, includ, list)
676 720 : DEALLOCATE (lrot, lwght)
677 :
678 720 : 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 720 : 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 720 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
695 :
696 720 : iou = csym%punit
697 720 : IF (iou > 0) THEN
698 : WRITE (iou, '(/,(T2,A79))') &
699 302 : "*******************************************************************************", &
700 302 : "** Special K-Point Generation by SPGLIB **", &
701 604 : "*******************************************************************************"
702 : END IF
703 720 : CALL cite_reference(Togo2018)
704 :
705 2160 : ALLOCATE (srot(3, 3, csym%n_operations))
706 720 : CALL setup_spglib_operations(csym, srot, nrot)
707 720 : CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
708 720 : DEALLOCATE (srot)
709 :
710 720 : IF (iou > 0) THEN
711 : WRITE (iou, '((T2,A79))') &
712 302 : "*******************************************************************************", &
713 302 : "** Finished SPGLIB **", &
714 604 : "*******************************************************************************"
715 : END IF
716 :
717 720 : 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 900 : 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 900 : 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 900 : CPASSERT(csym%symlib)
861 :
862 575448 : srot = 0
863 900 : csym%nrtot = 0
864 900 : IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
865 900 : IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
866 900 : IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
867 900 : IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
868 4500 : ALLOCATE (csym%rt(3, 3, csym%n_operations), csym%vt(3, csym%n_operations))
869 5400 : ALLOCATE (csym%ibrot(csym%n_operations), csym%f0(csym%nat, csym%n_operations))
870 575448 : csym%rt = 0.0_dp
871 177684 : csym%vt = 0.0_dp
872 45096 : csym%ibrot = 0
873 383840 : csym%f0 = 0
874 :
875 900 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
876 900 : h_inv = inv_3x3(csym%hmat)
877 2700 : ALLOCATE (perm(csym%nat))
878 :
879 900 : eye = 0
880 900 : eye(1, 1) = 1
881 900 : eye(2, 2) = 1
882 900 : eye(3, 3) = 1
883 :
884 900 : 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 4500 : DO pass = 1, 4
888 181284 : DO iop = 1, csym%n_operations
889 2298192 : irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
890 2298192 : frot(1:3, 1:3) = TRANSPOSE(irot(1:3, 1:3))
891 410840 : identity = ALL(frot == eye)
892 : zero_translation = ALL(ABS(csym%translations(1:3, iop) - &
893 283576 : ANINT(csym%translations(1:3, iop))) <= eps)
894 176784 : IF (pass == 1 .AND. (.NOT. identity .OR. .NOT. zero_translation)) CYCLE
895 133488 : IF (pass == 2 .AND. (identity .OR. .NOT. zero_translation)) CYCLE
896 95566 : IF (pass == 3 .AND. (.NOT. identity .OR. zero_translation)) CYCLE
897 52122 : IF (pass == 4 .AND. (identity .OR. zero_translation)) CYCLE
898 :
899 44196 : duplicate = .FALSE.
900 997734 : DO jop = 1, nrot
901 2488006 : IF (ALL(frot == srot(:, :, jop))) THEN
902 : duplicate = .TRUE.
903 : EXIT
904 : END IF
905 : END DO
906 44196 : IF (duplicate) CYCLE
907 :
908 13188 : CALL spglib_atom_permutation(csym, frot, csym%translations(:, iop), perm, valid)
909 13188 : IF (.NOT. valid) CYCLE
910 :
911 13188 : nrot = nrot + 1
912 :
913 171444 : srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
914 171444 : rfrac(1:3, 1:3) = REAL(frot(1:3, 1:3), KIND=dp)
915 1200108 : csym%rt(1:3, 1:3, nrot) = MATMUL(csym%hmat, MATMUL(rfrac, h_inv))
916 52752 : csym%vt(1:3, nrot) = csym%translations(1:3, iop)
917 13188 : csym%ibrot(nrot) = nrot
918 109772 : csym%f0(1:csym%nat, nrot) = perm(1:csym%nat)
919 : END DO
920 : END DO
921 :
922 900 : DEALLOCATE (perm)
923 900 : csym%nrtot = nrot
924 900 : IF (nrot == 0) CALL cp_abort(__LOCATION__, "SPGLIB did not return usable symmetry operations")
925 :
926 900 : END SUBROUTINE setup_spglib_operations
927 :
928 : ! **************************************************************************************************
929 : !> \brief Store unique SPGLIB rotations for K290-backend diagnostic reduction
930 : !> \param csym ...
931 : !> \param srot integer rotations in fractional coordinates
932 : !> \param nrot number of stored rotations
933 : ! **************************************************************************************************
934 20 : SUBROUTINE setup_spglib_reduction_rotations(csym, srot, nrot)
935 : TYPE(csym_type) :: csym
936 : INTEGER, DIMENSION(:, :, :), INTENT(OUT) :: srot
937 : INTEGER, INTENT(OUT) :: nrot
938 :
939 : INTEGER :: iop, jop, pass
940 : INTEGER, DIMENSION(3, 3) :: eye, frot, irot
941 : LOGICAL :: duplicate, identity
942 :
943 20 : CPASSERT(csym%symlib)
944 :
945 30388 : srot = 0
946 20 : eye = 0
947 20 : eye(1, 1) = 1
948 20 : eye(2, 2) = 1
949 20 : eye(3, 3) = 1
950 :
951 20 : nrot = 0
952 : ! Keep the identity first, matching the representative k-point operation.
953 60 : DO pass = 1, 2
954 4732 : DO iop = 1, csym%n_operations
955 60736 : irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
956 60736 : frot(1:3, 1:3) = TRANSPOSE(irot(1:3, 1:3))
957 10096 : identity = ALL(frot == eye)
958 4672 : IF (pass == 1 .AND. .NOT. identity) CYCLE
959 2392 : IF (pass == 2 .AND. identity) CYCLE
960 :
961 2336 : duplicate = .FALSE.
962 56528 : DO jop = 1, nrot
963 140872 : IF (ALL(frot == srot(:, :, jop))) THEN
964 : duplicate = .TRUE.
965 : EXIT
966 : END IF
967 : END DO
968 2336 : IF (duplicate) CYCLE
969 :
970 608 : nrot = nrot + 1
971 9728 : srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
972 : END DO
973 : END DO
974 :
975 20 : IF (nrot == 0) CALL cp_abort(__LOCATION__, "SPGLIB did not return usable symmetry rotations")
976 :
977 20 : END SUBROUTINE setup_spglib_reduction_rotations
978 :
979 : ! **************************************************************************************************
980 : !> \brief Determine the atom permutation generated by a SPGLIB space-group operation
981 : !> \param csym ...
982 : !> \param rot integer rotation in fractional coordinates
983 : !> \param trans fractional translation
984 : !> \param perm atom permutation
985 : !> \param valid whether all atoms were mapped
986 : ! **************************************************************************************************
987 13188 : SUBROUTINE spglib_atom_permutation(csym, rot, trans, perm, valid)
988 : TYPE(csym_type) :: csym
989 : INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
990 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: trans
991 : INTEGER, DIMENSION(:), INTENT(OUT) :: perm
992 : LOGICAL, INTENT(OUT) :: valid
993 :
994 : INTEGER :: i, j, nat
995 : LOGICAL :: found
996 13188 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
997 : REAL(KIND=dp) :: eps
998 : REAL(KIND=dp), DIMENSION(3) :: diff, spos
999 : REAL(KIND=dp), DIMENSION(3, 3) :: rfrac
1000 :
1001 13188 : nat = csym%nat
1002 13188 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1003 171444 : rfrac(1:3, 1:3) = REAL(rot(1:3, 1:3), KIND=dp)
1004 39564 : ALLOCATE (used(nat))
1005 13188 : used = .FALSE.
1006 106172 : perm = 0
1007 13188 : valid = .TRUE.
1008 :
1009 106172 : DO i = 1, nat
1010 1487744 : spos(1:3) = MATMUL(rfrac(1:3, 1:3), csym%scoord(1:3, i)) + trans(1:3)
1011 427076 : found = .FALSE.
1012 427076 : DO j = 1, nat
1013 427076 : IF (used(j)) CYCLE
1014 231706 : IF (csym%atype(i) /= csym%atype(j)) CYCLE
1015 926824 : diff(1:3) = spos(1:3) - csym%scoord(1:3, j)
1016 926824 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1017 528434 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1018 92984 : perm(i) = j
1019 92984 : used(j) = .TRUE.
1020 : found = .TRUE.
1021 : EXIT
1022 : END IF
1023 : END DO
1024 13188 : IF (.NOT. found) THEN
1025 0 : valid = .FALSE.
1026 0 : EXIT
1027 : END IF
1028 : END DO
1029 :
1030 13188 : DEALLOCATE (used)
1031 :
1032 13188 : END SUBROUTINE spglib_atom_permutation
1033 :
1034 : ! **************************************************************************************************
1035 : !> \brief Reduce a k-point mesh with SPGLIB direct-space operations
1036 : !> \param csym ...
1037 : !> \param xkp full k-point mesh in reciprocal lattice coordinates
1038 : !> \param wkp reduced k-point weights
1039 : !> \param kpop symmetry operation mapping the representative k-point to a mesh point
1040 : !> \param srot integer rotations in fractional coordinates
1041 : !> \param nrot number of stored rotations
1042 : ! **************************************************************************************************
1043 882 : SUBROUTINE reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
1044 : TYPE(csym_type) :: csym
1045 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1046 : REAL(KIND=dp), DIMENSION(:) :: wkp
1047 : INTEGER, DIMENSION(:) :: kpop
1048 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
1049 : INTEGER, INTENT(IN) :: nrot
1050 :
1051 : INTEGER :: i, iop, isign, j, kr, nkpts, score
1052 882 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kscore
1053 : INTEGER, DIMENSION(3, 3) :: krot
1054 : REAL(KIND=dp) :: eps
1055 : REAL(KIND=dp), DIMENSION(3) :: diff, rr
1056 :
1057 882 : nkpts = SIZE(wkp)
1058 882 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1059 2646 : ALLOCATE (kscore(nkpts))
1060 :
1061 13110 : wkp = 0.0_dp
1062 13110 : kpop = 0
1063 13110 : csym%kplink(1, :) = 0
1064 13110 : kscore = HUGE(0)
1065 :
1066 13110 : DO i = 1, nkpts
1067 12228 : IF (csym%kplink(1, i) /= 0) CYCLE
1068 :
1069 2122 : csym%kplink(1, i) = i
1070 2122 : wkp(i) = 1.0_dp
1071 2122 : kpop(i) = 1
1072 2122 : kscore(i) = 0
1073 :
1074 35606 : DO iop = 1, nrot
1075 32602 : kr = csym%ibrot(iop)
1076 32602 : krot = reciprocal_rotation(srot(:, :, kr))
1077 32602 : score = spglib_operation_score(csym, iop, srot(:, :, kr))
1078 110034 : DO isign = 1, 2
1079 1630100 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp(1:3, i))
1080 65204 : IF (isign == 2) THEN
1081 130408 : rr(1:3) = -rr(1:3)
1082 32602 : kr = -csym%ibrot(iop)
1083 : ELSE
1084 32602 : kr = csym%ibrot(iop)
1085 : END IF
1086 :
1087 1764606 : DO j = 1, nkpts
1088 7047096 : diff(1:3) = xkp(1:3, j) - rr(1:3)
1089 7047096 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1090 2414970 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1091 62372 : IF (csym%kplink(1, j) == 0) THEN
1092 10106 : csym%kplink(1, j) = i
1093 10106 : wkp(i) = wkp(i) + 1.0_dp
1094 10106 : kpop(j) = kr
1095 10106 : kscore(j) = score
1096 : ELSE
1097 52266 : CPASSERT(csym%kplink(1, j) == i)
1098 52266 : IF (score < kscore(j)) THEN
1099 3304 : kpop(j) = kr
1100 3304 : kscore(j) = score
1101 : END IF
1102 : END IF
1103 : EXIT
1104 : END IF
1105 : END DO
1106 97806 : IF (j > nkpts) CYCLE
1107 : END DO
1108 : END DO
1109 : END DO
1110 :
1111 13110 : DO i = 1, nkpts
1112 12228 : CPASSERT(csym%kplink(1, i) /= 0)
1113 13110 : CPASSERT(kpop(i) /= 0)
1114 : END DO
1115 882 : DEALLOCATE (kscore)
1116 :
1117 882 : END SUBROUTINE reduce_spglib_kpoint_mesh
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief Reduce an explicit k-point set by inversion/time-reversal.
1121 : !> \param csym ...
1122 : !> \param xkp_full explicit k-point coordinates
1123 : !> \param wkp_full explicit k-point weights
1124 : ! **************************************************************************************************
1125 0 : SUBROUTINE reduce_general_inversion(csym, xkp_full, wkp_full)
1126 : TYPE(csym_type) :: csym
1127 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1128 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_full
1129 :
1130 : INTEGER :: i, j, nfull, nred
1131 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1132 0 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1133 : REAL(KIND=dp) :: eps
1134 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1135 : REAL(KIND=dp), DIMENSION(3) :: diff
1136 :
1137 0 : nfull = SIZE(wkp_full)
1138 0 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1139 0 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1140 0 : used = .FALSE.
1141 0 : rep = 0
1142 0 : wred = 0.0_dp
1143 0 : nred = 0
1144 :
1145 0 : DO i = 1, nfull
1146 0 : IF (used(i)) CYCLE
1147 0 : nred = nred + 1
1148 0 : rep(nred) = i
1149 0 : used(i) = .TRUE.
1150 0 : csym%kplink(1, i) = i
1151 0 : csym%kpop(i) = 1
1152 0 : wred(nred) = wkp_full(i)
1153 0 : DO j = i + 1, nfull
1154 0 : IF (used(j)) CYCLE
1155 0 : diff(1:3) = xkp_full(1:3, j) + xkp_full(1:3, i)
1156 0 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1157 0 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1158 0 : IF (ABS(wkp_full(j) - wkp_full(i)) > eps) THEN
1159 : CALL cp_abort(__LOCATION__, &
1160 : "KPOINTS%INVERSION_SYMMETRY_ONLY with SCHEME GENERAL requires "// &
1161 0 : "equal weights for inversion-related k-points.")
1162 : END IF
1163 0 : used(j) = .TRUE.
1164 0 : csym%kplink(1, j) = i
1165 0 : csym%kpop(j) = -1
1166 0 : wred(nred) = wred(nred) + wkp_full(j)
1167 : END IF
1168 : END DO
1169 : END DO
1170 :
1171 0 : csym%nkpoint = nred
1172 0 : ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1173 0 : DO i = 1, nred
1174 0 : csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1175 0 : csym%wkpoint(i) = wred(i)
1176 : END DO
1177 0 : DO i = 1, nfull
1178 0 : DO j = 1, nred
1179 0 : IF (csym%kplink(1, i) == rep(j)) THEN
1180 0 : csym%kplink(2, i) = j
1181 0 : EXIT
1182 : END IF
1183 : END DO
1184 : END DO
1185 :
1186 0 : DEALLOCATE (rep, used, wred)
1187 :
1188 0 : END SUBROUTINE reduce_general_inversion
1189 :
1190 : ! **************************************************************************************************
1191 : !> \brief Reduce an explicit k-point set with K290 symmetry operations.
1192 : !> \param csym ...
1193 : !> \param xkp_full explicit k-point coordinates
1194 : ! **************************************************************************************************
1195 4 : SUBROUTINE reduce_general_k290(csym, xkp_full)
1196 : TYPE(csym_type) :: csym
1197 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1198 :
1199 : INTEGER :: i, ibsign, iop, j, kr, nfull, nred
1200 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1201 : LOGICAL :: found
1202 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1203 : REAL(KIND=dp) :: alat, eps
1204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1205 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, diff, rr, wcart
1206 :
1207 4 : nfull = SIZE(xkp_full, 2)
1208 4 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1209 4 : CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1210 :
1211 24 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1212 4 : used = .FALSE.
1213 4 : rep = 0
1214 4 : wred = 0.0_dp
1215 4 : nred = 0
1216 :
1217 36 : DO i = 1, nfull
1218 32 : IF (used(i)) CYCLE
1219 4 : nred = nred + 1
1220 4 : rep(nred) = i
1221 4 : used(i) = .TRUE.
1222 4 : csym%kplink(1, i) = i
1223 4 : csym%kpop(i) = 1
1224 4 : wred(nred) = 1.0_dp
1225 :
1226 200 : DO iop = 1, csym%nrtot
1227 608 : DO ibsign = 1, 2
1228 384 : kr = csym%ibrot(iop)
1229 : wcart(1:3) = alat*(xkp_full(1, i)*b1(1:3) + &
1230 : xkp_full(2, i)*b2(1:3) + &
1231 1536 : xkp_full(3, i)*b3(1:3))
1232 1536 : wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1233 384 : IF (ibsign == 2) THEN
1234 768 : wcart(1:3) = -wcart(1:3)
1235 192 : kr = -kr
1236 : END IF
1237 1536 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1238 1536 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1239 1536 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1240 :
1241 384 : found = .FALSE.
1242 1728 : DO j = 1, nfull
1243 6912 : diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1244 6912 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1245 3648 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1246 384 : found = .TRUE.
1247 384 : IF (.NOT. used(j)) THEN
1248 28 : used(j) = .TRUE.
1249 28 : csym%kplink(1, j) = i
1250 28 : csym%kpop(j) = kr
1251 28 : wred(nred) = wred(nred) + 1.0_dp
1252 : ELSE
1253 356 : CPASSERT(csym%kplink(1, j) == i)
1254 : END IF
1255 : EXIT
1256 : END IF
1257 : END DO
1258 192 : IF (.NOT. found) THEN
1259 : CALL cp_abort(__LOCATION__, &
1260 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1261 0 : "to be closed under the K290 symmetry operations.")
1262 : END IF
1263 : END DO
1264 : END DO
1265 : END DO
1266 :
1267 4 : CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1268 :
1269 4 : DEALLOCATE (rep, used, wred)
1270 :
1271 4 : END SUBROUTINE reduce_general_k290
1272 :
1273 : ! **************************************************************************************************
1274 : !> \brief Reduce an explicit k-point set with SPGLIB symmetry operations.
1275 : !> \param csym ...
1276 : !> \param xkp_full explicit k-point coordinates
1277 : ! **************************************************************************************************
1278 18 : SUBROUTINE reduce_general_spglib(csym, xkp_full)
1279 : TYPE(csym_type) :: csym
1280 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1281 :
1282 : INTEGER :: i, iop, isign, j, kr, nfull, nred, nrot
1283 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1284 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
1285 : INTEGER, DIMENSION(3, 3) :: krot
1286 : LOGICAL :: found
1287 18 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1288 : REAL(KIND=dp) :: eps
1289 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1290 : REAL(KIND=dp), DIMENSION(3) :: diff, rr
1291 :
1292 18 : nfull = SIZE(xkp_full, 2)
1293 18 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1294 54 : ALLOCATE (srot(3, 3, csym%n_operations))
1295 18 : CALL setup_spglib_operations(csym, srot, nrot)
1296 :
1297 108 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1298 18 : used = .FALSE.
1299 18 : rep = 0
1300 18 : wred = 0.0_dp
1301 18 : nred = 0
1302 :
1303 162 : DO i = 1, nfull
1304 144 : IF (used(i)) CYCLE
1305 18 : nred = nred + 1
1306 18 : rep(nred) = i
1307 18 : used(i) = .TRUE.
1308 18 : csym%kplink(1, i) = i
1309 18 : csym%kpop(i) = 1
1310 18 : wred(nred) = 1.0_dp
1311 :
1312 900 : DO iop = 1, nrot
1313 864 : kr = csym%ibrot(iop)
1314 864 : krot = reciprocal_rotation(srot(:, :, kr))
1315 2736 : DO isign = 1, 2
1316 43200 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp_full(1:3, i))
1317 1728 : IF (isign == 2) THEN
1318 3456 : rr(1:3) = -rr(1:3)
1319 864 : kr = -csym%ibrot(iop)
1320 : ELSE
1321 864 : kr = csym%ibrot(iop)
1322 : END IF
1323 :
1324 1728 : found = .FALSE.
1325 7776 : DO j = 1, nfull
1326 31104 : diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1327 31104 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1328 16416 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1329 1728 : found = .TRUE.
1330 1728 : IF (.NOT. used(j)) THEN
1331 126 : used(j) = .TRUE.
1332 126 : csym%kplink(1, j) = i
1333 126 : csym%kpop(j) = kr
1334 126 : wred(nred) = wred(nred) + 1.0_dp
1335 : ELSE
1336 1602 : CPASSERT(csym%kplink(1, j) == i)
1337 : END IF
1338 : EXIT
1339 : END IF
1340 : END DO
1341 864 : IF (.NOT. found) THEN
1342 : CALL cp_abort(__LOCATION__, &
1343 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1344 0 : "to be closed under the requested symmetry operations.")
1345 : END IF
1346 : END DO
1347 : END DO
1348 : END DO
1349 :
1350 18 : CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1351 :
1352 18 : DEALLOCATE (rep, srot, used, wred)
1353 :
1354 18 : END SUBROUTINE reduce_general_spglib
1355 :
1356 : ! **************************************************************************************************
1357 : !> \brief Reduce an explicit k-point set with SPGLIB rotations and K290 operations.
1358 : !> \param csym ...
1359 : !> \param xkp_full explicit k-point coordinates
1360 : ! **************************************************************************************************
1361 4 : SUBROUTINE reduce_general_spglib_k290(csym, xkp_full)
1362 : TYPE(csym_type) :: csym
1363 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1364 :
1365 : INTEGER :: i, iop, isign, j, k290_op, nfull, nred, &
1366 : nrot, nskipped
1367 : INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1368 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
1369 : INTEGER, DIMENSION(3, 3) :: krot
1370 : LOGICAL :: found, valid
1371 4 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1372 : REAL(KIND=dp) :: alat, eps
1373 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wred
1374 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, diff, rr
1375 :
1376 4 : nfull = SIZE(xkp_full, 2)
1377 4 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1378 4 : CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1379 12 : ALLOCATE (srot(3, 3, csym%n_operations))
1380 4 : CALL setup_spglib_reduction_rotations(csym, srot, nrot)
1381 :
1382 24 : ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1383 4 : used = .FALSE.
1384 4 : rep = 0
1385 4 : wred = 0.0_dp
1386 4 : nred = 0
1387 4 : nskipped = 0
1388 :
1389 36 : DO i = 1, nfull
1390 32 : IF (used(i)) CYCLE
1391 4 : nred = nred + 1
1392 4 : rep(nred) = i
1393 4 : used(i) = .TRUE.
1394 4 : csym%kplink(1, i) = i
1395 4 : csym%kpop(i) = 1
1396 4 : wred(nred) = 1.0_dp
1397 :
1398 200 : DO iop = 1, nrot
1399 192 : krot = reciprocal_rotation(srot(:, :, iop))
1400 608 : DO isign = 1, 2
1401 9600 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp_full(1:3, i))
1402 960 : IF (isign == 2) rr(1:3) = -rr(1:3)
1403 :
1404 1728 : found = .FALSE.
1405 1728 : DO j = 1, nfull
1406 6912 : diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1407 6912 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1408 3648 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1409 384 : found = .TRUE.
1410 : CALL find_k290_kpoint_operation(csym, xkp_full(1:3, i), xkp_full(1:3, j), &
1411 : a1, a2, a3, b1, b2, b3, alat, &
1412 384 : k290_op, valid)
1413 384 : IF (.NOT. valid) THEN
1414 0 : nskipped = nskipped + 1
1415 : EXIT
1416 : END IF
1417 384 : IF (.NOT. used(j)) THEN
1418 28 : used(j) = .TRUE.
1419 28 : csym%kplink(1, j) = i
1420 28 : csym%kpop(j) = k290_op
1421 28 : wred(nred) = wred(nred) + 1.0_dp
1422 : ELSE
1423 356 : CPASSERT(csym%kplink(1, j) == i)
1424 : END IF
1425 : EXIT
1426 : END IF
1427 : END DO
1428 192 : IF (.NOT. found) THEN
1429 : CALL cp_abort(__LOCATION__, &
1430 : "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1431 0 : "to be closed under the SPGLIB symmetry operations.")
1432 : END IF
1433 : END DO
1434 : END DO
1435 : END DO
1436 :
1437 4 : IF (nskipped > 0) THEN
1438 : CALL cp_warn(__LOCATION__, &
1439 : "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1440 0 : "the GENERAL k-point set was reduced only by the compatible mappings.")
1441 : END IF
1442 :
1443 4 : CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1444 :
1445 4 : DEALLOCATE (rep, srot, used, wred)
1446 :
1447 8 : END SUBROUTINE reduce_general_spglib_k290
1448 :
1449 : ! **************************************************************************************************
1450 : !> \brief Store reduced GENERAL k-point representatives.
1451 : !> \param csym ...
1452 : !> \param xkp_full explicit k-point coordinates
1453 : !> \param rep representative indices
1454 : !> \param wred representative multiplicities
1455 : !> \param nred number of reduced representatives
1456 : ! **************************************************************************************************
1457 26 : SUBROUTINE store_general_reduction(csym, xkp_full, rep, wred, nred)
1458 : TYPE(csym_type) :: csym
1459 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1460 : INTEGER, DIMENSION(:), INTENT(IN) :: rep
1461 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wred
1462 : INTEGER, INTENT(IN) :: nred
1463 :
1464 : INTEGER :: i, j, nfull
1465 :
1466 26 : nfull = SIZE(xkp_full, 2)
1467 26 : csym%nkpoint = nred
1468 130 : ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1469 52 : DO i = 1, nred
1470 104 : csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1471 52 : csym%wkpoint(i) = wred(i)
1472 : END DO
1473 234 : DO i = 1, nfull
1474 208 : DO j = 1, nred
1475 208 : IF (csym%kplink(1, i) == rep(j)) THEN
1476 208 : csym%kplink(2, i) = j
1477 208 : EXIT
1478 : END IF
1479 : END DO
1480 234 : CPASSERT(csym%kplink(2, i) /= 0)
1481 : END DO
1482 :
1483 26 : END SUBROUTINE store_general_reduction
1484 :
1485 : ! **************************************************************************************************
1486 : !> \brief Reduce a k-point mesh with SPGLIB rotations and K290 operations
1487 : !> \param csym ...
1488 : !> \param xkp full k-point mesh in reciprocal lattice coordinates
1489 : !> \param wkp reduced k-point weights
1490 : !> \param kpop K290 operation mapping the representative k-point to a mesh point
1491 : !> \param srot SPGLIB integer rotations in fractional coordinates
1492 : !> \param nrot number of stored rotations
1493 : !> \param a1 first lattice vector
1494 : !> \param a2 second lattice vector
1495 : !> \param a3 third lattice vector
1496 : !> \param b1 first reciprocal lattice vector
1497 : !> \param b2 second reciprocal lattice vector
1498 : !> \param b3 third reciprocal lattice vector
1499 : !> \param alat lattice scaling used by K290
1500 : ! **************************************************************************************************
1501 16 : SUBROUTINE reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
1502 : a1, a2, a3, b1, b2, b3, alat)
1503 : TYPE(csym_type) :: csym
1504 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1505 : REAL(KIND=dp), DIMENSION(:) :: wkp
1506 : INTEGER, DIMENSION(:) :: kpop
1507 : INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
1508 : INTEGER, INTENT(IN) :: nrot
1509 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1510 : REAL(KIND=dp), INTENT(IN) :: alat
1511 :
1512 : INTEGER :: i, iop, isign, j, k290_op, nkpts, &
1513 : nskipped
1514 : INTEGER, DIMENSION(3, 3) :: krot
1515 : LOGICAL :: valid
1516 : REAL(KIND=dp) :: eps
1517 : REAL(KIND=dp), DIMENSION(3) :: diff, rr
1518 :
1519 16 : nkpts = SIZE(wkp)
1520 16 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1521 16 : nskipped = 0
1522 :
1523 368 : wkp = 0.0_dp
1524 368 : kpop = 0
1525 368 : csym%kplink(1, :) = 0
1526 :
1527 368 : DO i = 1, nkpts
1528 352 : IF (csym%kplink(1, i) /= 0) CYCLE
1529 :
1530 64 : csym%kplink(1, i) = i
1531 64 : wkp(i) = 1.0_dp
1532 64 : kpop(i) = 1
1533 :
1534 688 : DO iop = 1, nrot
1535 608 : krot = reciprocal_rotation(srot(:, :, iop))
1536 2176 : DO isign = 1, 2
1537 30400 : rr(1:3) = MATMUL(REAL(krot(1:3, 1:3), KIND=dp), xkp(1:3, i))
1538 3040 : IF (isign == 2) rr(1:3) = -rr(1:3)
1539 :
1540 16224 : DO j = 1, nkpts
1541 64896 : diff(1:3) = xkp(1:3, j) - rr(1:3)
1542 64896 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1543 25184 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1544 1216 : IF (j == i) EXIT
1545 1024 : IF (csym%kplink(1, j) /= 0) THEN
1546 736 : CPASSERT(csym%kplink(1, j) == i)
1547 : EXIT
1548 : END IF
1549 :
1550 : CALL find_k290_kpoint_operation(csym, xkp(1:3, i), xkp(1:3, j), &
1551 : a1, a2, a3, b1, b2, b3, alat, &
1552 288 : k290_op, valid)
1553 288 : IF (.NOT. valid) THEN
1554 0 : nskipped = nskipped + 1
1555 0 : EXIT
1556 : END IF
1557 288 : csym%kplink(1, j) = i
1558 288 : wkp(i) = wkp(i) + 1.0_dp
1559 288 : kpop(j) = k290_op
1560 288 : EXIT
1561 : END IF
1562 : END DO
1563 1824 : IF (j > nkpts) CYCLE
1564 : END DO
1565 : END DO
1566 : END DO
1567 :
1568 368 : DO i = 1, nkpts
1569 352 : CPASSERT(csym%kplink(1, i) /= 0)
1570 368 : CPASSERT(kpop(i) /= 0)
1571 : END DO
1572 16 : IF (nskipped > 0) THEN
1573 : CALL cp_warn(__LOCATION__, &
1574 : "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1575 0 : "the mesh was reduced only by the compatible mappings.")
1576 : END IF
1577 :
1578 16 : END SUBROUTINE reduce_spglib_kpoint_mesh_k290
1579 :
1580 : ! **************************************************************************************************
1581 : !> \brief Find a K290 operation that maps one fractional k-point to another
1582 : !> \param csym ...
1583 : !> \param xref representative k-point
1584 : !> \param xtarget target k-point
1585 : !> \param a1 first lattice vector
1586 : !> \param a2 second lattice vector
1587 : !> \param a3 third lattice vector
1588 : !> \param b1 first reciprocal lattice vector
1589 : !> \param b2 second reciprocal lattice vector
1590 : !> \param b3 third reciprocal lattice vector
1591 : !> \param alat lattice scaling used by K290
1592 : !> \param k290_op K290 operation identifier
1593 : !> \param valid whether a matching K290 operation was found
1594 : ! **************************************************************************************************
1595 672 : SUBROUTINE find_k290_kpoint_operation(csym, xref, xtarget, a1, a2, a3, b1, b2, b3, alat, &
1596 : k290_op, valid)
1597 : TYPE(csym_type) :: csym
1598 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xref, xtarget, a1, a2, a3, b1, b2, b3
1599 : REAL(KIND=dp), INTENT(IN) :: alat
1600 : INTEGER, INTENT(OUT) :: k290_op
1601 : LOGICAL, INTENT(OUT) :: valid
1602 :
1603 : INTEGER :: ibsign, iop, kr
1604 : REAL(KIND=dp) :: eps
1605 : REAL(KIND=dp), DIMENSION(3) :: diff, rr, wcart
1606 :
1607 672 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1608 672 : k290_op = 0
1609 672 : valid = .FALSE.
1610 :
1611 1616 : DO iop = 1, csym%nrtot
1612 1616 : IF (iop > SIZE(csym%rt, 3)) CYCLE
1613 1616 : IF (csym%ibrot(iop) == 0) CYCLE
1614 3872 : DO ibsign = 1, 2
1615 11712 : wcart(1:3) = alat*(xref(1)*b1(1:3) + xref(2)*b2(1:3) + xref(3)*b3(1:3))
1616 11712 : wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1617 2928 : IF (ibsign == 2) THEN
1618 5248 : wcart(1:3) = -wcart(1:3)
1619 1312 : kr = -csym%ibrot(iop)
1620 : ELSE
1621 1616 : kr = csym%ibrot(iop)
1622 : END IF
1623 11712 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1624 11712 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1625 11712 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1626 :
1627 11712 : diff(1:3) = xtarget(1:3) - rr(1:3)
1628 11712 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1629 7056 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1630 672 : k290_op = kr
1631 672 : valid = .TRUE.
1632 672 : RETURN
1633 : END IF
1634 : END DO
1635 : END DO
1636 :
1637 : END SUBROUTINE find_k290_kpoint_operation
1638 :
1639 : ! **************************************************************************************************
1640 : !> \brief Score SPGLIB operations to choose stable atom transformations
1641 : !> \param csym ...
1642 : !> \param iop operation index
1643 : !> \param srot integer rotation in fractional coordinates
1644 : !> \return score, lower values are preferred
1645 : ! **************************************************************************************************
1646 32602 : FUNCTION spglib_operation_score(csym, iop, srot) RESULT(score)
1647 : TYPE(csym_type), INTENT(IN) :: csym
1648 : INTEGER, INTENT(IN) :: iop
1649 : INTEGER, DIMENSION(3, 3), INTENT(IN) :: srot
1650 : INTEGER :: score
1651 :
1652 : INTEGER :: i, nat
1653 : INTEGER, DIMENSION(3, 3) :: eye, r2
1654 : REAL(KIND=dp) :: eps
1655 :
1656 32602 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1657 32602 : nat = SIZE(csym%f0, 1)
1658 32602 : score = 0
1659 280624 : DO i = 1, nat
1660 280624 : IF (csym%f0(i, iop) /= i) score = score + 100
1661 : END DO
1662 85240 : IF (ANY(ABS(csym%vt(1:3, iop) - ANINT(csym%vt(1:3, iop))) > eps)) score = score + 10
1663 :
1664 32602 : eye = 0
1665 32602 : eye(1, 1) = 1
1666 32602 : eye(2, 2) = 1
1667 32602 : eye(3, 3) = 1
1668 1304080 : r2(1:3, 1:3) = MATMUL(srot(1:3, 1:3), srot(1:3, 1:3))
1669 239054 : IF (ANY(r2(1:3, 1:3) /= eye(1:3, 1:3))) score = score + 1
1670 :
1671 32602 : END FUNCTION spglib_operation_score
1672 :
1673 : ! **************************************************************************************************
1674 : !> \brief Reciprocal-space rotation corresponding to a fractional direct-space rotation
1675 : !> \param rot direct-space rotation
1676 : !> \return reciprocal-space rotation
1677 : ! **************************************************************************************************
1678 34266 : FUNCTION reciprocal_rotation(rot) RESULT(krot)
1679 : INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
1680 : INTEGER, DIMENSION(3, 3) :: krot
1681 :
1682 : REAL(KIND=dp), DIMENSION(3, 3) :: rinv
1683 :
1684 445458 : rinv = inv_3x3(REAL(rot(1:3, 1:3), KIND=dp))
1685 445458 : krot(1:3, 1:3) = NINT(TRANSPOSE(rinv(1:3, 1:3)))
1686 :
1687 34266 : END FUNCTION reciprocal_rotation
1688 :
1689 : ! **************************************************************************************************
1690 : !> \brief Reduce a CP2K Monkhorst-Pack mesh using K290 symmetry operations
1691 : !> \param csym ...
1692 : !> \param xkp full k-point mesh in reciprocal lattice coordinates
1693 : !> \param wkp reduced k-point weights
1694 : !> \param kpop symmetry operation mapping the representative k-point to a mesh point
1695 : !> \param nc number of point group operations
1696 : !> \param ib K290 operation identifiers
1697 : !> \param r K290 rotation matrices
1698 : !> \param a1 first lattice vector
1699 : !> \param a2 second lattice vector
1700 : !> \param a3 third lattice vector
1701 : !> \param b1 first reciprocal lattice vector
1702 : !> \param b2 second reciprocal lattice vector
1703 : !> \param b3 third reciprocal lattice vector
1704 : !> \param alat lattice scaling used by K290
1705 : ! **************************************************************************************************
1706 542 : SUBROUTINE reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
1707 : TYPE(csym_type) :: csym
1708 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1709 : REAL(KIND=dp), DIMENSION(:) :: wkp
1710 : INTEGER, DIMENSION(:) :: kpop
1711 : INTEGER, INTENT(IN) :: nc
1712 : INTEGER, DIMENSION(48), INTENT(IN) :: ib
1713 : REAL(KIND=dp), DIMENSION(3, 3, 48), INTENT(IN) :: r
1714 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1715 : REAL(KIND=dp), INTENT(IN) :: alat
1716 :
1717 : INTEGER :: i, ibsign, iop, j, kr, nkpts, &
1718 : nskipped_overlaps
1719 : REAL(KIND=dp) :: eps
1720 : REAL(KIND=dp), DIMENSION(3) :: diff, rr, wcart
1721 :
1722 542 : nkpts = SIZE(wkp)
1723 542 : eps = MAX(1.e-12_dp, 10.0_dp*csym%delta)
1724 542 : nskipped_overlaps = 0
1725 :
1726 4048 : wkp = 0.0_dp
1727 4048 : kpop = 0
1728 4048 : csym%kplink(1, :) = 0
1729 :
1730 4048 : DO i = 1, nkpts
1731 3506 : IF (csym%kplink(1, i) /= 0) CYCLE
1732 :
1733 1156 : csym%kplink(1, i) = i
1734 1156 : wkp(i) = 1.0_dp
1735 1156 : kpop(i) = 1
1736 :
1737 5304 : DO iop = 1, nc
1738 14324 : DO ibsign = 1, 2
1739 7212 : kr = ib(iop)
1740 28848 : wcart(1:3) = alat*(xkp(1, i)*b1(1:3) + xkp(2, i)*b2(1:3) + xkp(3, i)*b3(1:3))
1741 28848 : wcart(1:3) = kp_apply_operation(wcart(1:3), r(1:3, 1:3, kr))
1742 7212 : IF (ibsign == 2) THEN
1743 14424 : wcart(1:3) = -wcart(1:3)
1744 3606 : kr = -kr
1745 : END IF
1746 28848 : rr(1) = DOT_PRODUCT(a1(1:3), wcart(1:3))/alat
1747 28848 : rr(2) = DOT_PRODUCT(a2(1:3), wcart(1:3))/alat
1748 28848 : rr(3) = DOT_PRODUCT(a3(1:3), wcart(1:3))/alat
1749 :
1750 62810 : DO j = 1, nkpts
1751 250360 : diff(1:3) = xkp(1:3, j) - rr(1:3)
1752 250360 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1753 104206 : IF (ALL(ABS(diff(1:3)) < eps)) THEN
1754 6992 : IF (csym%kplink(1, j) == 0) THEN
1755 2350 : csym%kplink(1, j) = i
1756 2350 : wkp(i) = wkp(i) + 1.0_dp
1757 2350 : kpop(j) = kr
1758 4642 : ELSE IF (csym%kplink(1, j) /= i) THEN
1759 : ! Approximate K290 operation sets need not be closed for structures whose
1760 : ! coordinates lie close to several symmetry tolerances. Keep the existing
1761 : ! disjoint orbit instead of aborting or double-counting this mesh point.
1762 8 : nskipped_overlaps = nskipped_overlaps + 1
1763 : END IF
1764 : EXIT
1765 : END IF
1766 : END DO
1767 : ! Some point-group operations are incompatible with the requested Monkhorst-Pack mesh.
1768 3606 : IF (j > nkpts) CYCLE
1769 : END DO
1770 : END DO
1771 : END DO
1772 :
1773 542 : IF (nskipped_overlaps > 0) THEN
1774 56 : wkp = 1.0_dp
1775 2 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
1776 56 : DO i = 1, nkpts
1777 56 : IF (wkp(i) > 0.0_dp) THEN
1778 28 : kpop(i) = 1
1779 : ELSE
1780 26 : kpop(i) = 2
1781 : END IF
1782 : END DO
1783 2 : csym%nrtot = 0
1784 2 : csym%inversion_only = .TRUE.
1785 : CALL cp_warn(__LOCATION__, &
1786 : "The K290 k-point operations produced overlapping, non-closed mesh orbits; "// &
1787 : "falling back to inversion/time-reversal symmetry. Use SYMMETRY_BACKEND SPGLIB "// &
1788 2 : "for full reduction by a closed crystallographic symmetry group.")
1789 : ELSE
1790 3992 : DO i = 1, nkpts
1791 3452 : CPASSERT(csym%kplink(1, i) /= 0)
1792 3992 : CPASSERT(kpop(i) /= 0)
1793 : END DO
1794 : END IF
1795 :
1796 542 : END SUBROUTINE reduce_kpoint_mesh
1797 :
1798 : !> \brief ...
1799 : !> \param nk ...
1800 : !> \param xkp ...
1801 : !> \param wkp ...
1802 : !> \param shift ...
1803 : !> \param gamma_centered ...
1804 : ! **************************************************************************************************
1805 3180 : SUBROUTINE full_grid_gen(nk, xkp, wkp, shift, gamma_centered)
1806 : INTEGER, INTENT(IN) :: nk(3)
1807 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
1808 : REAL(KIND=dp), DIMENSION(:) :: wkp
1809 : REAL(KIND=dp), INTENT(IN) :: shift(3)
1810 : LOGICAL, INTENT(IN), OPTIONAL :: gamma_centered
1811 :
1812 : INTEGER :: i, idim, ix, iy, iz
1813 : INTEGER, DIMENSION(3) :: ik
1814 : LOGICAL :: gamma_mesh
1815 : REAL(KIND=dp) :: kpt_latt(3)
1816 :
1817 3180 : IF (PRESENT(gamma_centered)) THEN
1818 3180 : gamma_mesh = gamma_centered
1819 : ELSE
1820 : gamma_mesh = .FALSE.
1821 : END IF
1822 :
1823 34198 : wkp = 0.0_dp
1824 3180 : i = 0
1825 9284 : DO ix = 1, nk(1)
1826 22214 : DO iy = 1, nk(2)
1827 50052 : DO iz = 1, nk(3)
1828 31018 : i = i + 1
1829 31018 : ik(1) = ix
1830 31018 : ik(2) = iy
1831 31018 : ik(3) = iz
1832 124072 : DO idim = 1, 3
1833 124072 : IF (gamma_mesh .AND. MOD(nk(idim), 2) == 0) THEN
1834 : kpt_latt(idim) = REAL(2*ik(idim) - nk(idim), KIND=dp)/ &
1835 1216 : (2._dp*REAL(nk(idim), KIND=dp))
1836 : ELSE
1837 : kpt_latt(idim) = REAL(2*ik(idim) - nk(idim) - 1, KIND=dp)/ &
1838 91838 : (2._dp*REAL(nk(idim), KIND=dp))
1839 : END IF
1840 : END DO
1841 124072 : xkp(1:3, i) = kpt_latt(1:3)
1842 43948 : wkp(i) = 1.0_dp
1843 : END DO
1844 : END DO
1845 : END DO
1846 34198 : DO i = 1, nk(1)*nk(2)*nk(3)
1847 127252 : xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
1848 : END DO
1849 :
1850 3180 : END SUBROUTINE full_grid_gen
1851 :
1852 : ! **************************************************************************************************
1853 : !> \brief ...
1854 : !> \param xkp ...
1855 : !> \param wkp ...
1856 : !> \param link ...
1857 : ! **************************************************************************************************
1858 1672 : SUBROUTINE inversion_symm(xkp, wkp, link)
1859 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
1860 : REAL(KIND=dp), DIMENSION(:) :: wkp
1861 : INTEGER, DIMENSION(:) :: link
1862 :
1863 : INTEGER :: i, j, nkpts
1864 : REAL(KIND=dp), DIMENSION(3) :: diff
1865 :
1866 1672 : nkpts = SIZE(wkp, 1)
1867 :
1868 15858 : link(:) = 0
1869 15858 : DO i = 1, nkpts
1870 14186 : IF (link(i) == 0) link(i) = i
1871 173546 : DO j = i + 1, nkpts
1872 164556 : IF (wkp(j) == 0) CYCLE
1873 449440 : diff(1:3) = xkp(1:3, i) + xkp(1:3, j)
1874 449440 : diff(1:3) = diff(1:3) - ANINT(diff(1:3))
1875 169972 : IF (ALL(ABS(diff(1:3)) < 1.e-12_dp)) THEN
1876 6868 : wkp(i) = wkp(i) + wkp(j)
1877 6868 : wkp(j) = 0.0_dp
1878 6868 : link(j) = i
1879 6868 : EXIT
1880 : END IF
1881 : END DO
1882 : END DO
1883 :
1884 1672 : END SUBROUTINE inversion_symm
1885 :
1886 : ! **************************************************************************************************
1887 : !> \brief ...
1888 : !> \param x ...
1889 : !> \param r ...
1890 : !> \return ...
1891 : ! **************************************************************************************************
1892 10524 : FUNCTION kp_apply_operation(x, r) RESULT(y)
1893 : REAL(KIND=dp), INTENT(IN) :: x(3), r(3, 3)
1894 : REAL(KIND=dp) :: y(3)
1895 :
1896 10524 : y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
1897 10524 : y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
1898 10524 : y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
1899 :
1900 10524 : END FUNCTION kp_apply_operation
1901 :
1902 : ! **************************************************************************************************
1903 : !> \brief ...
1904 : !> \param csym ...
1905 : ! **************************************************************************************************
1906 3319 : SUBROUTINE print_crys_symmetry(csym)
1907 : TYPE(csym_type) :: csym
1908 :
1909 : INTEGER :: i, iunit, j, plevel
1910 :
1911 3319 : iunit = csym%punit
1912 3319 : IF (iunit >= 0) THEN
1913 1464 : plevel = csym%plevel
1914 1464 : WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
1915 1464 : IF (csym%symlib) THEN
1916 1462 : WRITE (iunit, '(A,T71,A10)') " International Symbol: ", ADJUSTR(TRIM(csym%international_symbol))
1917 1462 : WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", ADJUSTR(TRIM(csym%pointgroup_symbol))
1918 1462 : WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", ADJUSTR(TRIM(csym%schoenflies))
1919 : !
1920 1462 : WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
1921 1462 : IF (plevel > 0) THEN
1922 0 : DO i = 1, csym%n_operations
1923 : WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
1924 0 : " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
1925 0 : WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
1926 : END DO
1927 : END IF
1928 : ELSE
1929 2 : IF (csym%spglib_requested) THEN
1930 1 : WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not available"
1931 : ELSE
1932 1 : WRITE (iunit, "(T2,A)") "SPGLIB Crystal Symmetry Information was not requested"
1933 : END IF
1934 : END IF
1935 : END IF
1936 :
1937 3319 : END SUBROUTINE print_crys_symmetry
1938 :
1939 : ! **************************************************************************************************
1940 : !> \brief ...
1941 : !> \param csym ...
1942 : ! **************************************************************************************************
1943 3034 : SUBROUTINE print_kp_symmetry(csym)
1944 : TYPE(csym_type), INTENT(IN) :: csym
1945 :
1946 : INTEGER :: i, iunit, nat, nmesh, plevel
1947 :
1948 3034 : iunit = csym%punit
1949 3034 : IF (iunit >= 0) THEN
1950 1179 : plevel = csym%plevel
1951 1179 : WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
1952 1179 : WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
1953 1179 : WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
1954 4618 : DO i = 1, csym%nkpoint
1955 4618 : WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), NINT(csym%wkpoint(i))
1956 : END DO
1957 1179 : nmesh = csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
1958 1179 : IF (nmesh > 0) THEN
1959 1171 : WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
1960 : ELSE
1961 8 : nmesh = SIZE(csym%kpmesh, 2)
1962 8 : WRITE (iunit, '(/,A,T70,I10)') " Explicit K-point Set: ", nmesh
1963 : END IF
1964 1179 : WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
1965 10805 : DO i = 1, nmesh
1966 9626 : WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
1967 20431 : csym%kplink(1:2, i), csym%kpop(i)
1968 : END DO
1969 1179 : IF (csym%nrtot > 0) THEN
1970 600 : WRITE (iunit, '(/,A)') " Atom Transformation Table"
1971 600 : nat = SIZE(csym%f0, 1)
1972 5507 : DO i = 1, csym%nrtot
1973 5507 : WRITE (iunit, '(T10,A,I5,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
1974 : END DO
1975 : END IF
1976 : END IF
1977 :
1978 3034 : END SUBROUTINE print_kp_symmetry
1979 :
1980 0 : END MODULE cryssym
|
