(git:39fe4c4)
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 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
83CONTAINS
84
85! **************************************************************************************************
86!> \brief Release the CSYM type
87!> \param csym The CSYM type
88! **************************************************************************************************
89 SUBROUTINE release_csym_type(csym)
90 TYPE(csym_type) :: csym
91
92 IF (ALLOCATED(csym%rotations)) THEN
93 DEALLOCATE (csym%rotations)
94 END IF
95 IF (ALLOCATED(csym%translations)) THEN
96 DEALLOCATE (csym%translations)
97 END IF
98 IF (ALLOCATED(csym%atype)) THEN
99 DEALLOCATE (csym%atype)
100 END IF
101 IF (ALLOCATED(csym%scoord)) THEN
102 DEALLOCATE (csym%scoord)
103 END IF
104 IF (ALLOCATED(csym%xkpoint)) THEN
105 DEALLOCATE (csym%xkpoint)
106 END IF
107 IF (ALLOCATED(csym%wkpoint)) THEN
108 DEALLOCATE (csym%wkpoint)
109 END IF
110 IF (ALLOCATED(csym%kpmesh)) THEN
111 DEALLOCATE (csym%kpmesh)
112 END IF
113 IF (ALLOCATED(csym%kplink)) THEN
114 DEALLOCATE (csym%kplink)
115 END IF
116 IF (ALLOCATED(csym%kpop)) THEN
117 DEALLOCATE (csym%kpop)
118 END IF
119 IF (ALLOCATED(csym%rt)) THEN
120 DEALLOCATE (csym%rt)
121 END IF
122 IF (ALLOCATED(csym%vt)) THEN
123 DEALLOCATE (csym%vt)
124 END IF
125 IF (ALLOCATED(csym%f0)) THEN
126 DEALLOCATE (csym%f0)
127 END IF
128 IF (ALLOCATED(csym%ibrot)) THEN
129 DEALLOCATE (csym%ibrot)
130 END IF
131
132 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 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 CALL timeset(routinen, handle)
160
161 !..total number of atoms
162 nat = SIZE(scoor, 2)
163 csym%nat = nat
164
165 ! output unit
166 IF (PRESENT(iounit)) THEN
167 csym%punit = iounit
168 ELSE
169 csym%punit = -1
170 END IF
171
172 ! accuracy for symmetry
173 IF (PRESENT(delta)) THEN
174 csym%delta = delta
175 ELSE
176 csym%delta = 1.e-6_dp
177 END IF
178
179 !..set cell values
180 csym%hmat = hmat
181
182 ! atom types
183 ALLOCATE (csym%atype(nat))
184 csym%atype(1:nat) = types(1:nat)
185
186 ! scaled coordinates
187 ALLOCATE (csym%scoord(3, nat))
188 csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
189
190 csym%n_operations = 0
191
192 !..try spglib
193 my_use_spglib = .true.
194 IF (PRESENT(use_spglib)) my_use_spglib = use_spglib
195 csym%spglib_requested = my_use_spglib
196 IF (.NOT. my_use_spglib) THEN
197 spglib = .false.
198 ELSE
199 major = spg_get_major_version()
200 minor = spg_get_minor_version()
201 micro = spg_get_micro_version()
202 IF (major == 0) THEN
203 CALL cp_warn(__location__, "Symmetry library SPGLIB not available")
204 spglib = .false.
205 ELSE
206 spglib = .true.
207 CALL cite_reference(togo2018)
208 ierr = spg_get_international(csym%international_symbol, transpose(hmat), scoor, types, nat, delta)
209 IF (ierr == 0) THEN
210 CALL cp_warn(__location__, "Symmetry Library SPGLIB failed")
211 spglib = .false.
212 ELSE
213 nop = spg_get_multiplicity(transpose(hmat), scoor, types, nat, delta)
214 ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
215 csym%n_operations = nop
216 ierr = spg_get_symmetry(csym%rotations, csym%translations, nop, &
217 transpose(hmat), scoor, types, nat, delta)
218 ! Schoenflies Symbol
219 csym%schoenflies = ' '
220 ierr = spg_get_schoenflies(csym%schoenflies, transpose(hmat), scoor, types, nat, delta)
221 ! Point Group
222 csym%pointgroup_symbol = ' '
223 tra_mat = 0
224 ierr = spg_get_pointgroup(csym%pointgroup_symbol, tra_mat, &
225 csym%rotations, csym%n_operations)
226
227 CALL strip_control_codes(csym%international_symbol)
228 CALL strip_control_codes(csym%schoenflies)
229 CALL strip_control_codes(csym%pointgroup_symbol)
230 END IF
231 END IF
232 END IF
233 csym%symlib = spglib
234
235 CALL timestop(handle)
236
237 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 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 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 CALL timeset(routinen, handle)
272
273 IF (PRESENT(shift)) THEN
274 csym%wvk0 = shift
275 ELSE
276 csym%wvk0 = 0.0_dp
277 END IF
278
279 csym%istriz = -1
280 IF (PRESENT(symm)) THEN
281 IF (symm) csym%istriz = 1
282 END IF
283
284 IF (PRESENT(full_grid)) THEN
285 fullmesh = full_grid
286 ELSE
287 fullmesh = .false.
288 END IF
289 csym%fullgrid = fullmesh
290
291 IF (PRESENT(gamma_centered)) THEN
292 gamma_mesh = gamma_centered
293 ELSE
294 gamma_mesh = .false.
295 END IF
296
297 IF (PRESENT(inversion_symmetry_only)) THEN
298 inversion_only = inversion_symmetry_only
299 ELSE
300 inversion_only = .false.
301 END IF
302 csym%inversion_only = inversion_only
303
304 IF (PRESENT(use_spglib_reduction)) THEN
305 spglib_reduction = use_spglib_reduction
306 ELSE
307 spglib_reduction = .false.
308 END IF
309 csym%spglib_reduction = spglib_reduction
310
311 IF (PRESENT(use_spglib_backend)) THEN
312 spglib_backend = use_spglib_backend
313 ELSE
314 spglib_backend = .false.
315 END IF
316 csym%spglib_backend = spglib_backend
317
318 IF (spglib_backend .AND. .NOT. spglib_reduction) THEN
319 CALL cp_abort(__location__, &
320 "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 (spglib_backend .OR. spglib_reduction) .AND. .NOT. csym%symlib) THEN
324 CALL cp_abort(__location__, &
325 "SPGLIB k-point symmetry was requested, but SPGLIB is not available")
326 END IF
327
328 csym%nkpoint = 0
329 csym%mesh(1:3) = nk(1:3)
330 csym%nrtot = 0
331 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
332 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
333 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
334 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
335 ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
336
337 nkpts = nk(1)*nk(2)*nk(3)
338 ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
339 ! kp: link
340 ALLOCATE (csym%kplink(2, nkpts))
341 csym%kplink = 0
342 kpop = 0
343
344 ! go through all the options
345 IF (csym%symlib) THEN
346 ! symmetry library is available
347 IF (fullmesh) THEN
348 ! full mesh requested
349 CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
350 IF (csym%istriz == 1) THEN
351 ! use inversion symmetry
352 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
353 ELSE
354 ! full kpoint mesh is used
355 END IF
356 ELSE IF (csym%istriz /= 1 .OR. inversion_only) THEN
357 ! use inversion symmetry
358 CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
359 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
360 ELSE
361 ! use symmetry library to reduce k-points
362 CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
363 IF (spglib_backend) THEN
364 CALL kp_symmetry_spglib(csym, xkp, wkp, kpop)
365 ELSE
366 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 CALL full_grid_gen(nk, xkp, wkp, shift, gamma_centered=gamma_mesh)
373 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 CALL kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction=.false.)
376 ELSE IF (csym%istriz /= 1 .AND. fullmesh) THEN
377 ! full kpoint mesh is used
378 DO i = 1, nkpts
379 csym%kplink(1, i) = i
380 END DO
381 ELSE
382 ! use inversion symmetry
383 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
384 END IF
385 END IF
386 ! count kpoints
387 nkp = 0
388 DO i = 1, nkpts
389 IF (wkp(i) > 0.0_dp) nkp = nkp + 1
390 END DO
391
392 ! store reduced kpoint set
393 csym%nkpoint = nkp
394 ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
395 ALLOCATE (xptr(nkp))
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 xptr(j) = ik
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 DO j = 1, nkp
415 IF (i == xptr(j)) THEN
416 csym%kplink(2, ik) = j
417 EXIT
418 END IF
419 END DO
420 END DO
421 DEALLOCATE (xptr)
422
423 ! kp: operations
424 ALLOCATE (csym%kpop(nkpts))
425 IF (csym%nrtot > 0 .AND. csym%istriz == 1 .AND. .NOT. fullmesh .AND. &
426 .NOT. inversion_only) THEN
427 ! atomic symmetry operations possible
428 csym%kpop(1:nkpts) = kpop(1:nkpts)
429 DO ik = 1, nkpts
430 cpassert(csym%kpop(ik) /= 0)
431 END DO
432 ELSE
433 ! only time reversal symmetry
434 DO ik = 1, nkpts
435 IF (wkp(ik) > 0.0_dp) THEN
436 csym%kpop(ik) = 1
437 ELSE
438 csym%kpop(ik) = 2
439 END IF
440 END DO
441 END IF
442
443 DEALLOCATE (xkp, wkp, kpop)
444
445 CALL timestop(handle)
446
447 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 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 CALL timeset(routinen, handle)
479
480 nfull = SIZE(wkp_in)
481 cpassert(SIZE(xkp_in, 1) == 3)
482 cpassert(SIZE(xkp_in, 2) == nfull)
483
484 atomic_symmetry = .false.
485 IF (PRESENT(symm)) atomic_symmetry = symm
486 csym%istriz = -1
487 IF (atomic_symmetry) csym%istriz = 1
488 fullmesh = .false.
489 IF (PRESENT(full_grid)) fullmesh = full_grid
490 inversion_only = .false.
491 IF (PRESENT(inversion_symmetry_only)) inversion_only = inversion_symmetry_only
492 spglib_reduction = .false.
493 IF (PRESENT(use_spglib_reduction)) spglib_reduction = use_spglib_reduction
494 spglib_backend = .false.
495 IF (PRESENT(use_spglib_backend)) spglib_backend = use_spglib_backend
496
497 csym%fullgrid = fullmesh
498 csym%inversion_only = inversion_only
499 csym%spglib_reduction = spglib_reduction
500 csym%spglib_backend = spglib_backend
501 csym%nkpoint = 0
502 csym%mesh(1:3) = 0
503 csym%nrtot = 0
504 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
505 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
506 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
507 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
508 ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
509 IF (ALLOCATED(csym%xkpoint)) DEALLOCATE (csym%xkpoint)
510 IF (ALLOCATED(csym%wkpoint)) DEALLOCATE (csym%wkpoint)
511 IF (ALLOCATED(csym%kpmesh)) DEALLOCATE (csym%kpmesh)
512 IF (ALLOCATED(csym%kplink)) DEALLOCATE (csym%kplink)
513 IF (ALLOCATED(csym%kpop)) DEALLOCATE (csym%kpop)
514
515 ALLOCATE (csym%kpmesh(3, nfull), csym%kplink(2, nfull), csym%kpop(nfull))
516 csym%kpmesh(1:3, 1:nfull) = xkp_in(1:3, 1:nfull)
517 csym%kplink = 0
518 csym%kpop = 1
519
520 IF (.NOT. atomic_symmetry .OR. fullmesh) THEN
521 csym%nkpoint = nfull
522 ALLOCATE (csym%xkpoint(3, nfull), csym%wkpoint(nfull))
523 csym%xkpoint(1:3, 1:nfull) = xkp_in(1:3, 1:nfull)
524 csym%wkpoint(1:nfull) = wkp_in(1:nfull)
525 DO i = 1, nfull
526 csym%kplink(1:2, i) = i
527 END DO
528 ELSE IF (inversion_only) THEN
529 CALL reduce_general_inversion(csym, xkp_in, wkp_in)
530 ELSE
531 weight_eps = max(1.e-12_dp, 10.0_dp*csym%delta)
532 IF (any(abs(wkp_in(1:nfull) - wkp_in(1)) > weight_eps)) THEN
533 CALL cp_abort(__location__, &
534 "KPOINTS%SYMMETRY with SCHEME GENERAL requires equal explicit weights.")
535 END IF
536 IF (spglib_backend) THEN
537 IF (.NOT. csym%symlib) THEN
538 CALL cp_abort(__location__, &
539 "SCHEME GENERAL with SYMMETRY_BACKEND SPGLIB requires SPGLIB.")
540 END IF
541 CALL reduce_general_spglib(csym, xkp_in)
542 ELSE IF (spglib_reduction) THEN
543 IF (.NOT. csym%symlib) THEN
544 CALL cp_abort(__location__, &
545 "SCHEME GENERAL with SYMMETRY_REDUCTION_METHOD SPGLIB requires SPGLIB.")
546 END IF
547 CALL setup_k290_operations(csym)
548 CALL reduce_general_spglib_k290(csym, xkp_in)
549 ELSE
550 CALL setup_k290_operations(csym)
551 CALL reduce_general_k290(csym, xkp_in)
552 END IF
553 END IF
554
555 CALL timestop(handle)
556
557 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: includ, isc, list, lwght, ty
578 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0, lrot
579 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
580 INTEGER, DIMENSION(48) :: ib
581 LOGICAL :: spglib_reduction
582 REAL(kind=dp) :: alat
583 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 iou = csym%punit
590 hmat = csym%hmat
591 nat = csym%nat
592 iq1 = csym%mesh(1)
593 iq2 = csym%mesh(2)
594 iq3 = csym%mesh(3)
595 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 wvk0 = 0.0_dp
599 istriz = csym%istriz
600 IF (PRESENT(use_spglib_reduction)) THEN
601 spglib_reduction = use_spglib_reduction
602 ELSE
603 spglib_reduction = .false.
604 END IF
605 a1(1:3) = hmat(1:3, 1)
606 a2(1:3) = hmat(1:3, 2)
607 a3(1:3) = hmat(1:3, 3)
608 alat = sqrt(sum(a1**2))
609 strain = 0.0_dp
610 ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, 200), ty(nat), isc(nat), f0(49, nat))
611 ty(1:nat) = csym%atype(1:nat)
612 nsp = maxval(ty)
613 DO i = 1, nat
614 xkapa(1:3, i) = matmul(hmat, csym%scoord(1:3, i))
615 END DO
616 nhash = max(1000, nkpoint)
617 ALLOCATE (wvkl(3, nkpoint), rlist(3, nkpoint), includ(nkpoint), list(nhash + nkpoint))
618 ALLOCATE (lrot(48, nkpoint), lwght(nkpoint))
619
620 IF (iou > 0) THEN
621 WRITE (iou, '(/,(T2,A79))') &
622 "*******************************************************************************", &
623 "** Special K-Point Generation by K290 **", &
624 "*******************************************************************************"
625 END IF
626 CALL cite_reference(worlton1972)
627 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 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 vt, f0, r, tvec, origin, rx, isc, csym%delta)
637
638 IF (iou > 0) THEN
639 WRITE (iou, '((T2,A79))') &
640 "*******************************************************************************", &
641 "** Finished K290 **", &
642 "*******************************************************************************"
643 END IF
644
645 csym%nrtot = nc
646 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
647 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
648 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
649 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
650 ALLOCATE (csym%rt(3, 3, nc), csym%vt(3, nc), csym%ibrot(nc))
651 csym%vt(1:3, 1:nc) = vt(1:3, 1:nc)
652 ALLOCATE (csym%f0(nat, nc))
653 DO i = 1, nc
654 csym%rt(1:3, 1:3, i) = r(1:3, 1:3, ib(i))
655 csym%f0(1:nat, i) = f0(i, 1:nat)
656 END DO
657 csym%ibrot(1:nc) = ib(1:nc)
658
659 IF (csym%n_operations > nc .AND. .NOT. spglib_reduction) THEN
660 IF (ALLOCATED(srot)) DEALLOCATE (srot)
661 ALLOCATE (srot(3, 3, csym%n_operations))
662 CALL setup_spglib_operations(csym, srot, nrot)
663 CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
664 DEALLOCATE (srot)
665 ELSE IF (spglib_reduction) THEN
666 ALLOCATE (srot(3, 3, csym%n_operations))
667 CALL setup_spglib_reduction_rotations(csym, srot, nrot)
668 CALL reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
669 a1, a2, a3, b1, b2, b3, alat)
670 DEALLOCATE (srot)
671 ELSE
672 CALL reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
673 END IF
674 DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
675 DEALLOCATE (wvkl, rlist, includ, list)
676 DEALLOCATE (lrot, lwght)
677
678 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 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 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
695
696 iou = csym%punit
697 IF (iou > 0) THEN
698 WRITE (iou, '(/,(T2,A79))') &
699 "*******************************************************************************", &
700 "** Special K-Point Generation by SPGLIB **", &
701 "*******************************************************************************"
702 END IF
703 CALL cite_reference(togo2018)
704
705 ALLOCATE (srot(3, 3, csym%n_operations))
706 CALL setup_spglib_operations(csym, srot, nrot)
707 CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
708 DEALLOCATE (srot)
709
710 IF (iou > 0) THEN
711 WRITE (iou, '((T2,A79))') &
712 "*******************************************************************************", &
713 "** Finished SPGLIB **", &
714 "*******************************************************************************"
715 END IF
716
717 END SUBROUTINE kp_symmetry_spglib
718
719! **************************************************************************************************
720!> \brief Store K290 atomic symmetry operations without reducing a generated mesh.
721!> \param csym ...
722! **************************************************************************************************
723 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: includ, isc, list, lwght, ty
730 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0, lrot
731 INTEGER, DIMENSION(48) :: ib
732 REAL(kind=dp) :: alat
733 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 iou = csym%punit
740 nat = csym%nat
741 CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
742 iq1 = max(1, csym%mesh(1))
743 iq2 = max(1, csym%mesh(2))
744 iq3 = max(1, csym%mesh(3))
745 nkpoint = max(10, 10*iq1*iq2*iq3)
746 strain = 0.0_dp
747 wvk0 = 0.0_dp
748
749 ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, 200), ty(nat), isc(nat), f0(49, nat))
750 ty(1:nat) = csym%atype(1:nat)
751 nsp = maxval(ty)
752 DO i = 1, nat
753 xkapa(1:3, i) = matmul(csym%hmat, csym%scoord(1:3, i))
754 END DO
755 nhash = max(1000, nkpoint)
756 ALLOCATE (wvkl(3, nkpoint), rlist(3, nkpoint), includ(nkpoint), list(nhash + nkpoint))
757 ALLOCATE (lrot(48, nkpoint), lwght(nkpoint))
758
759 IF (iou > 0) THEN
760 WRITE (iou, '(/,(T2,A79))') &
761 "*******************************************************************************", &
762 "** Special K-Point Generation by K290 **", &
763 "*******************************************************************************"
764 END IF
765 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 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 vt, f0, r, tvec, origin, rx, isc, csym%delta)
775
776 IF (iou > 0) THEN
777 WRITE (iou, '((T2,A79))') &
778 "*******************************************************************************", &
779 "** Finished K290 **", &
780 "*******************************************************************************"
781 END IF
782
783 csym%nrtot = nc
784 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
785 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
786 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
787 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
788 ALLOCATE (csym%rt(3, 3, nc), csym%vt(3, nc), csym%ibrot(nc))
789 csym%vt(1:3, 1:nc) = vt(1:3, 1:nc)
790 ALLOCATE (csym%f0(nat, nc))
791 DO i = 1, nc
792 csym%rt(1:3, 1:3, i) = r(1:3, 1:3, ib(i))
793 csym%f0(1:nat, i) = f0(i, 1:nat)
794 END DO
795 csym%ibrot(1:nc) = ib(1:nc)
796
797 DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
798 DEALLOCATE (wvkl, rlist, includ, list)
799 DEALLOCATE (lrot, lwght)
800
801 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 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 a1(1:3) = csym%hmat(1:3, 1)
822 a2(1:3) = csym%hmat(1:3, 2)
823 a3(1:3) = csym%hmat(1:3, 3)
824 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 a2(3)*a3(2)*a1(1) - a3(3)*a1(2)*a2(1)
828 volum = abs(volum)
829 b1(1) = (a2(2)*a3(3) - a2(3)*a3(2))/volum
830 b1(2) = (a2(3)*a3(1) - a2(1)*a3(3))/volum
831 b1(3) = (a2(1)*a3(2) - a2(2)*a3(1))/volum
832 b2(1) = (a3(2)*a1(3) - a3(3)*a1(2))/volum
833 b2(2) = (a3(3)*a1(1) - a3(1)*a1(3))/volum
834 b2(3) = (a3(1)*a1(2) - a3(2)*a1(1))/volum
835 b3(1) = (a1(2)*a2(3) - a1(3)*a2(2))/volum
836 b3(2) = (a1(3)*a2(1) - a1(1)*a2(3))/volum
837 b3(3) = (a1(1)*a2(2) - a1(2)*a2(1))/volum
838
839 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 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 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 cpassert(csym%symlib)
861
862 srot = 0
863 csym%nrtot = 0
864 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
865 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
866 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
867 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
868 ALLOCATE (csym%rt(3, 3, csym%n_operations), csym%vt(3, csym%n_operations))
869 ALLOCATE (csym%ibrot(csym%n_operations), csym%f0(csym%nat, csym%n_operations))
870 csym%rt = 0.0_dp
871 csym%vt = 0.0_dp
872 csym%ibrot = 0
873 csym%f0 = 0
874
875 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
876 h_inv = inv_3x3(csym%hmat)
877 ALLOCATE (perm(csym%nat))
878
879 eye = 0
880 eye(1, 1) = 1
881 eye(2, 2) = 1
882 eye(3, 3) = 1
883
884 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 DO pass = 1, 4
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 zero_translation = all(abs(csym%translations(1:3, iop) - &
893 anint(csym%translations(1:3, iop))) <= eps)
894 IF (pass == 1 .AND. (.NOT. identity .OR. .NOT. zero_translation)) cycle
895 IF (pass == 2 .AND. (identity .OR. .NOT. zero_translation)) cycle
896 IF (pass == 3 .AND. (.NOT. identity .OR. zero_translation)) cycle
897 IF (pass == 4 .AND. (identity .OR. zero_translation)) cycle
898
899 duplicate = .false.
900 DO jop = 1, nrot
901 IF (all(frot == srot(:, :, jop))) THEN
902 duplicate = .true.
903 EXIT
904 END IF
905 END DO
906 IF (duplicate) cycle
907
908 CALL spglib_atom_permutation(csym, frot, csym%translations(:, iop), perm, valid)
909 IF (.NOT. valid) cycle
910
911 nrot = nrot + 1
912
913 srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
914 rfrac(1:3, 1:3) = real(frot(1:3, 1:3), kind=dp)
915 csym%rt(1:3, 1:3, nrot) = matmul(csym%hmat, matmul(rfrac, h_inv))
916 csym%vt(1:3, nrot) = csym%translations(1:3, iop)
917 csym%ibrot(nrot) = nrot
918 csym%f0(1:csym%nat, nrot) = perm(1:csym%nat)
919 END DO
920 END DO
921
922 DEALLOCATE (perm)
923 csym%nrtot = nrot
924 IF (nrot == 0) CALL cp_abort(__location__, "SPGLIB did not return usable symmetry operations")
925
926 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 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 cpassert(csym%symlib)
944
945 srot = 0
946 eye = 0
947 eye(1, 1) = 1
948 eye(2, 2) = 1
949 eye(3, 3) = 1
950
951 nrot = 0
952 ! Keep the identity first, matching the representative k-point operation.
953 DO pass = 1, 2
954 DO iop = 1, csym%n_operations
955 irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
956 frot(1:3, 1:3) = transpose(irot(1:3, 1:3))
957 identity = all(frot == eye)
958 IF (pass == 1 .AND. .NOT. identity) cycle
959 IF (pass == 2 .AND. identity) cycle
960
961 duplicate = .false.
962 DO jop = 1, nrot
963 IF (all(frot == srot(:, :, jop))) THEN
964 duplicate = .true.
965 EXIT
966 END IF
967 END DO
968 IF (duplicate) cycle
969
970 nrot = nrot + 1
971 srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
972 END DO
973 END DO
974
975 IF (nrot == 0) CALL cp_abort(__location__, "SPGLIB did not return usable symmetry rotations")
976
977 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 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 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 nat = csym%nat
1002 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1003 rfrac(1:3, 1:3) = real(rot(1:3, 1:3), kind=dp)
1004 ALLOCATE (used(nat))
1005 used = .false.
1006 perm = 0
1007 valid = .true.
1008
1009 DO i = 1, nat
1010 spos(1:3) = matmul(rfrac(1:3, 1:3), csym%scoord(1:3, i)) + trans(1:3)
1011 found = .false.
1012 DO j = 1, nat
1013 IF (used(j)) cycle
1014 IF (csym%atype(i) /= csym%atype(j)) cycle
1015 diff(1:3) = spos(1:3) - csym%scoord(1:3, j)
1016 diff(1:3) = diff(1:3) - anint(diff(1:3))
1017 IF (all(abs(diff(1:3)) < eps)) THEN
1018 perm(i) = j
1019 used(j) = .true.
1020 found = .true.
1021 EXIT
1022 END IF
1023 END DO
1024 IF (.NOT. found) THEN
1025 valid = .false.
1026 EXIT
1027 END IF
1028 END DO
1029
1030 DEALLOCATE (used)
1031
1032 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 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 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 nkpts = SIZE(wkp)
1058 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1059 ALLOCATE (kscore(nkpts))
1060
1061 wkp = 0.0_dp
1062 kpop = 0
1063 csym%kplink(1, :) = 0
1064 kscore = huge(0)
1065
1066 DO i = 1, nkpts
1067 IF (csym%kplink(1, i) /= 0) cycle
1068
1069 csym%kplink(1, i) = i
1070 wkp(i) = 1.0_dp
1071 kpop(i) = 1
1072 kscore(i) = 0
1073
1074 DO iop = 1, nrot
1075 kr = csym%ibrot(iop)
1076 krot = reciprocal_rotation(srot(:, :, kr))
1077 score = spglib_operation_score(csym, iop, srot(:, :, kr))
1078 DO isign = 1, 2
1079 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp(1:3, i))
1080 IF (isign == 2) THEN
1081 rr(1:3) = -rr(1:3)
1082 kr = -csym%ibrot(iop)
1083 ELSE
1084 kr = csym%ibrot(iop)
1085 END IF
1086
1087 DO j = 1, nkpts
1088 diff(1:3) = xkp(1:3, j) - rr(1:3)
1089 diff(1:3) = diff(1:3) - anint(diff(1:3))
1090 IF (all(abs(diff(1:3)) < eps)) THEN
1091 IF (csym%kplink(1, j) == 0) THEN
1092 csym%kplink(1, j) = i
1093 wkp(i) = wkp(i) + 1.0_dp
1094 kpop(j) = kr
1095 kscore(j) = score
1096 ELSE
1097 cpassert(csym%kplink(1, j) == i)
1098 IF (score < kscore(j)) THEN
1099 kpop(j) = kr
1100 kscore(j) = score
1101 END IF
1102 END IF
1103 EXIT
1104 END IF
1105 END DO
1106 IF (j > nkpts) cycle
1107 END DO
1108 END DO
1109 END DO
1110
1111 DO i = 1, nkpts
1112 cpassert(csym%kplink(1, i) /= 0)
1113 cpassert(kpop(i) /= 0)
1114 END DO
1115 DEALLOCATE (kscore)
1116
1117 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1132 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1133 REAL(kind=dp) :: eps
1134 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wred
1135 REAL(kind=dp), DIMENSION(3) :: diff
1136
1137 nfull = SIZE(wkp_full)
1138 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1139 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1140 used = .false.
1141 rep = 0
1142 wred = 0.0_dp
1143 nred = 0
1144
1145 DO i = 1, nfull
1146 IF (used(i)) cycle
1147 nred = nred + 1
1148 rep(nred) = i
1149 used(i) = .true.
1150 csym%kplink(1, i) = i
1151 csym%kpop(i) = 1
1152 wred(nred) = wkp_full(i)
1153 DO j = i + 1, nfull
1154 IF (used(j)) cycle
1155 diff(1:3) = xkp_full(1:3, j) + xkp_full(1:3, i)
1156 diff(1:3) = diff(1:3) - anint(diff(1:3))
1157 IF (all(abs(diff(1:3)) < eps)) THEN
1158 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 "equal weights for inversion-related k-points.")
1162 END IF
1163 used(j) = .true.
1164 csym%kplink(1, j) = i
1165 csym%kpop(j) = -1
1166 wred(nred) = wred(nred) + wkp_full(j)
1167 END IF
1168 END DO
1169 END DO
1170
1171 csym%nkpoint = nred
1172 ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1173 DO i = 1, nred
1174 csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1175 csym%wkpoint(i) = wred(i)
1176 END DO
1177 DO i = 1, nfull
1178 DO j = 1, nred
1179 IF (csym%kplink(1, i) == rep(j)) THEN
1180 csym%kplink(2, i) = j
1181 EXIT
1182 END IF
1183 END DO
1184 END DO
1185
1186 DEALLOCATE (rep, used, wred)
1187
1188 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 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 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 nfull = SIZE(xkp_full, 2)
1208 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1209 CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1210
1211 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1212 used = .false.
1213 rep = 0
1214 wred = 0.0_dp
1215 nred = 0
1216
1217 DO i = 1, nfull
1218 IF (used(i)) cycle
1219 nred = nred + 1
1220 rep(nred) = i
1221 used(i) = .true.
1222 csym%kplink(1, i) = i
1223 csym%kpop(i) = 1
1224 wred(nred) = 1.0_dp
1225
1226 DO iop = 1, csym%nrtot
1227 DO ibsign = 1, 2
1228 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 xkp_full(3, i)*b3(1:3))
1232 wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1233 IF (ibsign == 2) THEN
1234 wcart(1:3) = -wcart(1:3)
1235 kr = -kr
1236 END IF
1237 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1238 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1239 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1240
1241 found = .false.
1242 DO j = 1, nfull
1243 diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1244 diff(1:3) = diff(1:3) - anint(diff(1:3))
1245 IF (all(abs(diff(1:3)) < eps)) THEN
1246 found = .true.
1247 IF (.NOT. used(j)) THEN
1248 used(j) = .true.
1249 csym%kplink(1, j) = i
1250 csym%kpop(j) = kr
1251 wred(nred) = wred(nred) + 1.0_dp
1252 ELSE
1253 cpassert(csym%kplink(1, j) == i)
1254 END IF
1255 EXIT
1256 END IF
1257 END DO
1258 IF (.NOT. found) THEN
1259 CALL cp_abort(__location__, &
1260 "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1261 "to be closed under the K290 symmetry operations.")
1262 END IF
1263 END DO
1264 END DO
1265 END DO
1266
1267 CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1268
1269 DEALLOCATE (rep, used, wred)
1270
1271 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 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 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 nfull = SIZE(xkp_full, 2)
1293 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1294 ALLOCATE (srot(3, 3, csym%n_operations))
1295 CALL setup_spglib_operations(csym, srot, nrot)
1296
1297 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1298 used = .false.
1299 rep = 0
1300 wred = 0.0_dp
1301 nred = 0
1302
1303 DO i = 1, nfull
1304 IF (used(i)) cycle
1305 nred = nred + 1
1306 rep(nred) = i
1307 used(i) = .true.
1308 csym%kplink(1, i) = i
1309 csym%kpop(i) = 1
1310 wred(nred) = 1.0_dp
1311
1312 DO iop = 1, nrot
1313 kr = csym%ibrot(iop)
1314 krot = reciprocal_rotation(srot(:, :, kr))
1315 DO isign = 1, 2
1316 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp_full(1:3, i))
1317 IF (isign == 2) THEN
1318 rr(1:3) = -rr(1:3)
1319 kr = -csym%ibrot(iop)
1320 ELSE
1321 kr = csym%ibrot(iop)
1322 END IF
1323
1324 found = .false.
1325 DO j = 1, nfull
1326 diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1327 diff(1:3) = diff(1:3) - anint(diff(1:3))
1328 IF (all(abs(diff(1:3)) < eps)) THEN
1329 found = .true.
1330 IF (.NOT. used(j)) THEN
1331 used(j) = .true.
1332 csym%kplink(1, j) = i
1333 csym%kpop(j) = kr
1334 wred(nred) = wred(nred) + 1.0_dp
1335 ELSE
1336 cpassert(csym%kplink(1, j) == i)
1337 END IF
1338 EXIT
1339 END IF
1340 END DO
1341 IF (.NOT. found) THEN
1342 CALL cp_abort(__location__, &
1343 "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1344 "to be closed under the requested symmetry operations.")
1345 END IF
1346 END DO
1347 END DO
1348 END DO
1349
1350 CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1351
1352 DEALLOCATE (rep, srot, used, wred)
1353
1354 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 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 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 nfull = SIZE(xkp_full, 2)
1377 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1378 CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1379 ALLOCATE (srot(3, 3, csym%n_operations))
1380 CALL setup_spglib_reduction_rotations(csym, srot, nrot)
1381
1382 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1383 used = .false.
1384 rep = 0
1385 wred = 0.0_dp
1386 nred = 0
1387 nskipped = 0
1388
1389 DO i = 1, nfull
1390 IF (used(i)) cycle
1391 nred = nred + 1
1392 rep(nred) = i
1393 used(i) = .true.
1394 csym%kplink(1, i) = i
1395 csym%kpop(i) = 1
1396 wred(nred) = 1.0_dp
1397
1398 DO iop = 1, nrot
1399 krot = reciprocal_rotation(srot(:, :, iop))
1400 DO isign = 1, 2
1401 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp_full(1:3, i))
1402 IF (isign == 2) rr(1:3) = -rr(1:3)
1403
1404 found = .false.
1405 DO j = 1, nfull
1406 diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1407 diff(1:3) = diff(1:3) - anint(diff(1:3))
1408 IF (all(abs(diff(1:3)) < eps)) THEN
1409 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 k290_op, valid)
1413 IF (.NOT. valid) THEN
1414 nskipped = nskipped + 1
1415 EXIT
1416 END IF
1417 IF (.NOT. used(j)) THEN
1418 used(j) = .true.
1419 csym%kplink(1, j) = i
1420 csym%kpop(j) = k290_op
1421 wred(nred) = wred(nred) + 1.0_dp
1422 ELSE
1423 cpassert(csym%kplink(1, j) == i)
1424 END IF
1425 EXIT
1426 END IF
1427 END DO
1428 IF (.NOT. found) THEN
1429 CALL cp_abort(__location__, &
1430 "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1431 "to be closed under the SPGLIB symmetry operations.")
1432 END IF
1433 END DO
1434 END DO
1435 END DO
1436
1437 IF (nskipped > 0) THEN
1438 CALL cp_warn(__location__, &
1439 "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1440 "the GENERAL k-point set was reduced only by the compatible mappings.")
1441 END IF
1442
1443 CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1444
1445 DEALLOCATE (rep, srot, used, wred)
1446
1447 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 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 nfull = SIZE(xkp_full, 2)
1467 csym%nkpoint = nred
1468 ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1469 DO i = 1, nred
1470 csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1471 csym%wkpoint(i) = wred(i)
1472 END DO
1473 DO i = 1, nfull
1474 DO j = 1, nred
1475 IF (csym%kplink(1, i) == rep(j)) THEN
1476 csym%kplink(2, i) = j
1477 EXIT
1478 END IF
1479 END DO
1480 cpassert(csym%kplink(2, i) /= 0)
1481 END DO
1482
1483 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 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 nkpts = SIZE(wkp)
1520 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1521 nskipped = 0
1522
1523 wkp = 0.0_dp
1524 kpop = 0
1525 csym%kplink(1, :) = 0
1526
1527 DO i = 1, nkpts
1528 IF (csym%kplink(1, i) /= 0) cycle
1529
1530 csym%kplink(1, i) = i
1531 wkp(i) = 1.0_dp
1532 kpop(i) = 1
1533
1534 DO iop = 1, nrot
1535 krot = reciprocal_rotation(srot(:, :, iop))
1536 DO isign = 1, 2
1537 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp(1:3, i))
1538 IF (isign == 2) rr(1:3) = -rr(1:3)
1539
1540 DO j = 1, nkpts
1541 diff(1:3) = xkp(1:3, j) - rr(1:3)
1542 diff(1:3) = diff(1:3) - anint(diff(1:3))
1543 IF (all(abs(diff(1:3)) < eps)) THEN
1544 IF (j == i) EXIT
1545 IF (csym%kplink(1, j) /= 0) THEN
1546 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 k290_op, valid)
1553 IF (.NOT. valid) THEN
1554 nskipped = nskipped + 1
1555 EXIT
1556 END IF
1557 csym%kplink(1, j) = i
1558 wkp(i) = wkp(i) + 1.0_dp
1559 kpop(j) = k290_op
1560 EXIT
1561 END IF
1562 END DO
1563 IF (j > nkpts) cycle
1564 END DO
1565 END DO
1566 END DO
1567
1568 DO i = 1, nkpts
1569 cpassert(csym%kplink(1, i) /= 0)
1570 cpassert(kpop(i) /= 0)
1571 END DO
1572 IF (nskipped > 0) THEN
1573 CALL cp_warn(__location__, &
1574 "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1575 "the mesh was reduced only by the compatible mappings.")
1576 END IF
1577
1578 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 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 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1608 k290_op = 0
1609 valid = .false.
1610
1611 DO iop = 1, csym%nrtot
1612 IF (iop > SIZE(csym%rt, 3)) cycle
1613 IF (csym%ibrot(iop) == 0) cycle
1614 DO ibsign = 1, 2
1615 wcart(1:3) = alat*(xref(1)*b1(1:3) + xref(2)*b2(1:3) + xref(3)*b3(1:3))
1616 wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1617 IF (ibsign == 2) THEN
1618 wcart(1:3) = -wcart(1:3)
1619 kr = -csym%ibrot(iop)
1620 ELSE
1621 kr = csym%ibrot(iop)
1622 END IF
1623 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1624 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1625 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1626
1627 diff(1:3) = xtarget(1:3) - rr(1:3)
1628 diff(1:3) = diff(1:3) - anint(diff(1:3))
1629 IF (all(abs(diff(1:3)) < eps)) THEN
1630 k290_op = kr
1631 valid = .true.
1632 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 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 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1657 nat = SIZE(csym%f0, 1)
1658 score = 0
1659 DO i = 1, nat
1660 IF (csym%f0(i, iop) /= i) score = score + 100
1661 END DO
1662 IF (any(abs(csym%vt(1:3, iop) - anint(csym%vt(1:3, iop))) > eps)) score = score + 10
1663
1664 eye = 0
1665 eye(1, 1) = 1
1666 eye(2, 2) = 1
1667 eye(3, 3) = 1
1668 r2(1:3, 1:3) = matmul(srot(1:3, 1:3), srot(1:3, 1:3))
1669 IF (any(r2(1:3, 1:3) /= eye(1:3, 1:3))) score = score + 1
1670
1671 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 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 rinv = inv_3x3(real(rot(1:3, 1:3), kind=dp))
1685 krot(1:3, 1:3) = nint(transpose(rinv(1:3, 1:3)))
1686
1687 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 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 nkpts = SIZE(wkp)
1723 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1724 nskipped_overlaps = 0
1725
1726 wkp = 0.0_dp
1727 kpop = 0
1728 csym%kplink(1, :) = 0
1729
1730 DO i = 1, nkpts
1731 IF (csym%kplink(1, i) /= 0) cycle
1732
1733 csym%kplink(1, i) = i
1734 wkp(i) = 1.0_dp
1735 kpop(i) = 1
1736
1737 DO iop = 1, nc
1738 DO ibsign = 1, 2
1739 kr = ib(iop)
1740 wcart(1:3) = alat*(xkp(1, i)*b1(1:3) + xkp(2, i)*b2(1:3) + xkp(3, i)*b3(1:3))
1741 wcart(1:3) = kp_apply_operation(wcart(1:3), r(1:3, 1:3, kr))
1742 IF (ibsign == 2) THEN
1743 wcart(1:3) = -wcart(1:3)
1744 kr = -kr
1745 END IF
1746 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1747 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1748 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1749
1750 DO j = 1, nkpts
1751 diff(1:3) = xkp(1:3, j) - rr(1:3)
1752 diff(1:3) = diff(1:3) - anint(diff(1:3))
1753 IF (all(abs(diff(1:3)) < eps)) THEN
1754 IF (csym%kplink(1, j) == 0) THEN
1755 csym%kplink(1, j) = i
1756 wkp(i) = wkp(i) + 1.0_dp
1757 kpop(j) = kr
1758 ELSEIF (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 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 IF (j > nkpts) cycle
1769 END DO
1770 END DO
1771 END DO
1772
1773 IF (nskipped_overlaps > 0) THEN
1774 wkp = 1.0_dp
1775 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
1776 DO i = 1, nkpts
1777 IF (wkp(i) > 0.0_dp) THEN
1778 kpop(i) = 1
1779 ELSE
1780 kpop(i) = 2
1781 END IF
1782 END DO
1783 csym%nrtot = 0
1784 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 "for full reduction by a closed crystallographic symmetry group.")
1789 ELSE
1790 DO i = 1, nkpts
1791 cpassert(csym%kplink(1, i) /= 0)
1792 cpassert(kpop(i) /= 0)
1793 END DO
1794 END IF
1795
1796 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 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 IF (PRESENT(gamma_centered)) THEN
1818 gamma_mesh = gamma_centered
1819 ELSE
1820 gamma_mesh = .false.
1821 END IF
1822
1823 wkp = 0.0_dp
1824 i = 0
1825 DO ix = 1, nk(1)
1826 DO iy = 1, nk(2)
1827 DO iz = 1, nk(3)
1828 i = i + 1
1829 ik(1) = ix
1830 ik(2) = iy
1831 ik(3) = iz
1832 DO idim = 1, 3
1833 IF (gamma_mesh .AND. mod(nk(idim), 2) == 0) THEN
1834 kpt_latt(idim) = real(2*ik(idim) - nk(idim), kind=dp)/ &
1835 (2._dp*real(nk(idim), kind=dp))
1836 ELSE
1837 kpt_latt(idim) = real(2*ik(idim) - nk(idim) - 1, kind=dp)/ &
1838 (2._dp*real(nk(idim), kind=dp))
1839 END IF
1840 END DO
1841 xkp(1:3, i) = kpt_latt(1:3)
1842 wkp(i) = 1.0_dp
1843 END DO
1844 END DO
1845 END DO
1846 DO i = 1, nk(1)*nk(2)*nk(3)
1847 xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
1848 END DO
1849
1850 END SUBROUTINE full_grid_gen
1851
1852! **************************************************************************************************
1853!> \brief ...
1854!> \param xkp ...
1855!> \param wkp ...
1856!> \param link ...
1857! **************************************************************************************************
1858 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 nkpts = SIZE(wkp, 1)
1867
1868 link(:) = 0
1869 DO i = 1, nkpts
1870 IF (link(i) == 0) link(i) = i
1871 DO j = i + 1, nkpts
1872 IF (wkp(j) == 0) cycle
1873 diff(1:3) = xkp(1:3, i) + xkp(1:3, j)
1874 diff(1:3) = diff(1:3) - anint(diff(1:3))
1875 IF (all(abs(diff(1:3)) < 1.e-12_dp)) THEN
1876 wkp(i) = wkp(i) + wkp(j)
1877 wkp(j) = 0.0_dp
1878 link(j) = i
1879 EXIT
1880 END IF
1881 END DO
1882 END DO
1883
1884 END SUBROUTINE inversion_symm
1885
1886! **************************************************************************************************
1887!> \brief ...
1888!> \param x ...
1889!> \param r ...
1890!> \return ...
1891! **************************************************************************************************
1892 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 y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
1897 y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
1898 y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
1899
1900 END FUNCTION kp_apply_operation
1901
1902! **************************************************************************************************
1903!> \brief ...
1904!> \param csym ...
1905! **************************************************************************************************
1906 SUBROUTINE print_crys_symmetry(csym)
1907 TYPE(csym_type) :: csym
1908
1909 INTEGER :: i, iunit, j, plevel
1910
1911 iunit = csym%punit
1912 IF (iunit >= 0) THEN
1913 plevel = csym%plevel
1914 WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
1915 IF (csym%symlib) THEN
1916 WRITE (iunit, '(A,T71,A10)') " International Symbol: ", adjustr(trim(csym%international_symbol))
1917 WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", adjustr(trim(csym%pointgroup_symbol))
1918 WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", adjustr(trim(csym%schoenflies))
1919 !
1920 WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
1921 IF (plevel > 0) THEN
1922 DO i = 1, csym%n_operations
1923 WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
1924 " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
1925 WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
1926 END DO
1927 END IF
1928 ELSE
1929 IF (csym%spglib_requested) THEN
1930 WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not available"
1931 ELSE
1932 WRITE (iunit, "(T2,A)") "SPGLIB Crystal Symmetry Information was not requested"
1933 END IF
1934 END IF
1935 END IF
1936
1937 END SUBROUTINE print_crys_symmetry
1938
1939! **************************************************************************************************
1940!> \brief ...
1941!> \param csym ...
1942! **************************************************************************************************
1943 SUBROUTINE print_kp_symmetry(csym)
1944 TYPE(csym_type), INTENT(IN) :: csym
1945
1946 INTEGER :: i, iunit, nat, nmesh, plevel
1947
1948 iunit = csym%punit
1949 IF (iunit >= 0) THEN
1950 plevel = csym%plevel
1951 WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
1952 WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
1953 WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
1954 DO i = 1, csym%nkpoint
1955 WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), nint(csym%wkpoint(i))
1956 END DO
1957 nmesh = csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
1958 IF (nmesh > 0) THEN
1959 WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
1960 ELSE
1961 nmesh = SIZE(csym%kpmesh, 2)
1962 WRITE (iunit, '(/,A,T70,I10)') " Explicit K-point Set: ", nmesh
1963 END IF
1964 WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
1965 DO i = 1, nmesh
1966 WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
1967 csym%kplink(1:2, i), csym%kpop(i)
1968 END DO
1969 IF (csym%nrtot > 0) THEN
1970 WRITE (iunit, '(/,A)') " Atom Transformation Table"
1971 nat = SIZE(csym%f0, 1)
1972 DO i = 1, csym%nrtot
1973 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 END SUBROUTINE print_kp_symmetry
1979
1980END MODULE cryssym
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:1907
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:253
cryssym::release_csym_type
subroutine, public release_csym_type(csym)
Release the CSYM type.
Definition cryssym.F:90
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:462
cryssym::print_kp_symmetry
subroutine, public print_kp_symmetry(csym)
...
Definition cryssym.F:1944
cryssym::crys_sym_gen
subroutine, public crys_sym_gen(csym, scoor, types, hmat, delta, iounit, use_spglib)
...
Definition cryssym.F:145
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpsym
K-points and crystal symmetry routines based on.
Definition kpsym.F:28
kpsym::k290s
subroutine, public k290s(iout, nat, nkpoint, nsp, iq1, iq2, iq3, istriz, a1, a2, a3, alat, strain, xkapa, rx, tvec, ty, isc, f0, ntvec, wvk0, wvkl, lwght, lrot, nhash, includ, list, rlist, delta)
...
Definition kpsym.F:82
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:561
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
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:44