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