(git:618e178)
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 = hmat(1, 1)
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 = 1000
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 = 1000
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 = csym%hmat(1, 1)
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)) .AND. &
902 all(abs(csym%translations(1:3, iop) - csym%vt(1:3, jop) - &
903 anint(csym%translations(1:3, iop) - csym%vt(1:3, jop))) <= eps)) THEN
904 duplicate = .true.
905 EXIT
906 END IF
907 END DO
908 IF (duplicate) cycle
909
910 CALL spglib_atom_permutation(csym, frot, csym%translations(:, iop), perm, valid)
911 IF (.NOT. valid) cycle
912
913 nrot = nrot + 1
914
915 srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
916 rfrac(1:3, 1:3) = real(frot(1:3, 1:3), kind=dp)
917 csym%rt(1:3, 1:3, nrot) = matmul(csym%hmat, matmul(rfrac, h_inv))
918 csym%vt(1:3, nrot) = csym%translations(1:3, iop)
919 csym%ibrot(nrot) = nrot
920 csym%f0(1:csym%nat, nrot) = perm(1:csym%nat)
921 END DO
922 END DO
923
924 DEALLOCATE (perm)
925 csym%nrtot = nrot
926 IF (nrot == 0) CALL cp_abort(__location__, "SPGLIB did not return usable symmetry operations")
927
928 END SUBROUTINE setup_spglib_operations
929
930! **************************************************************************************************
931!> \brief Store unique SPGLIB rotations for K290-backend diagnostic reduction
932!> \param csym ...
933!> \param srot integer rotations in fractional coordinates
934!> \param nrot number of stored rotations
935! **************************************************************************************************
936 SUBROUTINE setup_spglib_reduction_rotations(csym, srot, nrot)
937 TYPE(csym_type) :: csym
938 INTEGER, DIMENSION(:, :, :), INTENT(OUT) :: srot
939 INTEGER, INTENT(OUT) :: nrot
940
941 INTEGER :: iop, jop, pass
942 INTEGER, DIMENSION(3, 3) :: eye, frot, irot
943 LOGICAL :: duplicate, identity
944
945 cpassert(csym%symlib)
946
947 srot = 0
948 eye = 0
949 eye(1, 1) = 1
950 eye(2, 2) = 1
951 eye(3, 3) = 1
952
953 nrot = 0
954 ! Keep the identity first, matching the representative k-point operation.
955 DO pass = 1, 2
956 DO iop = 1, csym%n_operations
957 irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
958 frot(1:3, 1:3) = transpose(irot(1:3, 1:3))
959 identity = all(frot == eye)
960 IF (pass == 1 .AND. .NOT. identity) cycle
961 IF (pass == 2 .AND. identity) cycle
962
963 duplicate = .false.
964 DO jop = 1, nrot
965 IF (all(frot == srot(:, :, jop))) THEN
966 duplicate = .true.
967 EXIT
968 END IF
969 END DO
970 IF (duplicate) cycle
971
972 nrot = nrot + 1
973 srot(1:3, 1:3, nrot) = frot(1:3, 1:3)
974 END DO
975 END DO
976
977 IF (nrot == 0) CALL cp_abort(__location__, "SPGLIB did not return usable symmetry rotations")
978
979 END SUBROUTINE setup_spglib_reduction_rotations
980
981! **************************************************************************************************
982!> \brief Determine the atom permutation generated by a SPGLIB space-group operation
983!> \param csym ...
984!> \param rot integer rotation in fractional coordinates
985!> \param trans fractional translation
986!> \param perm atom permutation
987!> \param valid whether all atoms were mapped
988! **************************************************************************************************
989 SUBROUTINE spglib_atom_permutation(csym, rot, trans, perm, valid)
990 TYPE(csym_type) :: csym
991 INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
992 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: trans
993 INTEGER, DIMENSION(:), INTENT(OUT) :: perm
994 LOGICAL, INTENT(OUT) :: valid
995
996 INTEGER :: i, j, nat
997 LOGICAL :: found
998 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
999 REAL(kind=dp) :: eps
1000 REAL(kind=dp), DIMENSION(3) :: diff, spos
1001 REAL(kind=dp), DIMENSION(3, 3) :: rfrac
1002
1003 nat = csym%nat
1004 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1005 rfrac(1:3, 1:3) = real(rot(1:3, 1:3), kind=dp)
1006 ALLOCATE (used(nat))
1007 used = .false.
1008 perm = 0
1009 valid = .true.
1010
1011 DO i = 1, nat
1012 spos(1:3) = matmul(rfrac(1:3, 1:3), csym%scoord(1:3, i)) + trans(1:3)
1013 found = .false.
1014 DO j = 1, nat
1015 IF (used(j)) cycle
1016 IF (csym%atype(i) /= csym%atype(j)) cycle
1017 diff(1:3) = spos(1:3) - csym%scoord(1:3, j)
1018 diff(1:3) = diff(1:3) - anint(diff(1:3))
1019 IF (all(abs(diff(1:3)) < eps)) THEN
1020 perm(i) = j
1021 used(j) = .true.
1022 found = .true.
1023 EXIT
1024 END IF
1025 END DO
1026 IF (.NOT. found) THEN
1027 valid = .false.
1028 EXIT
1029 END IF
1030 END DO
1031
1032 DEALLOCATE (used)
1033
1034 END SUBROUTINE spglib_atom_permutation
1035
1036! **************************************************************************************************
1037!> \brief Reduce a k-point mesh with SPGLIB direct-space operations
1038!> \param csym ...
1039!> \param xkp full k-point mesh in reciprocal lattice coordinates
1040!> \param wkp reduced k-point weights
1041!> \param kpop symmetry operation mapping the representative k-point to a mesh point
1042!> \param srot integer rotations in fractional coordinates
1043!> \param nrot number of stored rotations
1044! **************************************************************************************************
1045 SUBROUTINE reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
1046 TYPE(csym_type) :: csym
1047 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1048 REAL(kind=dp), DIMENSION(:) :: wkp
1049 INTEGER, DIMENSION(:) :: kpop
1050 INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
1051 INTEGER, INTENT(IN) :: nrot
1052
1053 INTEGER :: i, iop, isign, j, kr, nkpts, score
1054 INTEGER, ALLOCATABLE, DIMENSION(:) :: kscore
1055 INTEGER, DIMENSION(3, 3) :: krot
1056 REAL(kind=dp) :: eps
1057 REAL(kind=dp), DIMENSION(3) :: diff, rr
1058
1059 nkpts = SIZE(wkp)
1060 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1061 ALLOCATE (kscore(nkpts))
1062
1063 wkp = 0.0_dp
1064 kpop = 0
1065 csym%kplink(1, :) = 0
1066 kscore = huge(0)
1067
1068 DO i = 1, nkpts
1069 IF (csym%kplink(1, i) /= 0) cycle
1070
1071 csym%kplink(1, i) = i
1072 wkp(i) = 1.0_dp
1073 kpop(i) = 1
1074 kscore(i) = 0
1075
1076 DO iop = 1, nrot
1077 kr = csym%ibrot(iop)
1078 krot = reciprocal_rotation(srot(:, :, kr))
1079 score = spglib_operation_score(csym, iop, srot(:, :, kr))
1080 DO isign = 1, 2
1081 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp(1:3, i))
1082 IF (isign == 2) THEN
1083 rr(1:3) = -rr(1:3)
1084 kr = -csym%ibrot(iop)
1085 ELSE
1086 kr = csym%ibrot(iop)
1087 END IF
1088
1089 DO j = 1, nkpts
1090 diff(1:3) = xkp(1:3, j) - rr(1:3)
1091 diff(1:3) = diff(1:3) - anint(diff(1:3))
1092 IF (all(abs(diff(1:3)) < eps)) THEN
1093 IF (csym%kplink(1, j) == 0) THEN
1094 csym%kplink(1, j) = i
1095 wkp(i) = wkp(i) + 1.0_dp
1096 kpop(j) = kr
1097 kscore(j) = score
1098 ELSE
1099 cpassert(csym%kplink(1, j) == i)
1100 IF (score < kscore(j)) THEN
1101 kpop(j) = kr
1102 kscore(j) = score
1103 END IF
1104 END IF
1105 EXIT
1106 END IF
1107 END DO
1108 IF (j > nkpts) cycle
1109 END DO
1110 END DO
1111 END DO
1112
1113 DO i = 1, nkpts
1114 cpassert(csym%kplink(1, i) /= 0)
1115 cpassert(kpop(i) /= 0)
1116 END DO
1117 DEALLOCATE (kscore)
1118
1119 END SUBROUTINE reduce_spglib_kpoint_mesh
1120
1121! **************************************************************************************************
1122!> \brief Reduce an explicit k-point set by inversion/time-reversal.
1123!> \param csym ...
1124!> \param xkp_full explicit k-point coordinates
1125!> \param wkp_full explicit k-point weights
1126! **************************************************************************************************
1127 SUBROUTINE reduce_general_inversion(csym, xkp_full, wkp_full)
1128 TYPE(csym_type) :: csym
1129 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1130 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wkp_full
1131
1132 INTEGER :: i, j, nfull, nred
1133 INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1134 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1135 REAL(kind=dp) :: eps
1136 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wred
1137 REAL(kind=dp), DIMENSION(3) :: diff
1138
1139 nfull = SIZE(wkp_full)
1140 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
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) = wkp_full(i)
1155 DO j = i + 1, nfull
1156 IF (used(j)) cycle
1157 diff(1:3) = xkp_full(1:3, j) + xkp_full(1:3, i)
1158 diff(1:3) = diff(1:3) - anint(diff(1:3))
1159 IF (all(abs(diff(1:3)) < eps)) THEN
1160 IF (abs(wkp_full(j) - wkp_full(i)) > eps) THEN
1161 CALL cp_abort(__location__, &
1162 "KPOINTS%INVERSION_SYMMETRY_ONLY with SCHEME GENERAL requires "// &
1163 "equal weights for inversion-related k-points.")
1164 END IF
1165 used(j) = .true.
1166 csym%kplink(1, j) = i
1167 csym%kpop(j) = -1
1168 wred(nred) = wred(nred) + wkp_full(j)
1169 END IF
1170 END DO
1171 END DO
1172
1173 csym%nkpoint = nred
1174 ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1175 DO i = 1, nred
1176 csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1177 csym%wkpoint(i) = wred(i)
1178 END DO
1179 DO i = 1, nfull
1180 DO j = 1, nred
1181 IF (csym%kplink(1, i) == rep(j)) THEN
1182 csym%kplink(2, i) = j
1183 EXIT
1184 END IF
1185 END DO
1186 END DO
1187
1188 DEALLOCATE (rep, used, wred)
1189
1190 END SUBROUTINE reduce_general_inversion
1191
1192! **************************************************************************************************
1193!> \brief Reduce an explicit k-point set with K290 symmetry operations.
1194!> \param csym ...
1195!> \param xkp_full explicit k-point coordinates
1196! **************************************************************************************************
1197 SUBROUTINE reduce_general_k290(csym, xkp_full)
1198 TYPE(csym_type) :: csym
1199 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1200
1201 INTEGER :: i, ibsign, iop, j, kr, nfull, nred
1202 INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1203 LOGICAL :: found
1204 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1205 REAL(kind=dp) :: alat, eps
1206 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wred
1207 REAL(kind=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, diff, rr, wcart
1208
1209 nfull = SIZE(xkp_full, 2)
1210 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1211 CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1212
1213 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1214 used = .false.
1215 rep = 0
1216 wred = 0.0_dp
1217 nred = 0
1218
1219 DO i = 1, nfull
1220 IF (used(i)) cycle
1221 nred = nred + 1
1222 rep(nred) = i
1223 used(i) = .true.
1224 csym%kplink(1, i) = i
1225 csym%kpop(i) = 1
1226 wred(nred) = 1.0_dp
1227
1228 DO iop = 1, csym%nrtot
1229 DO ibsign = 1, 2
1230 kr = csym%ibrot(iop)
1231 wcart(1:3) = alat*(xkp_full(1, i)*b1(1:3) + &
1232 xkp_full(2, i)*b2(1:3) + &
1233 xkp_full(3, i)*b3(1:3))
1234 wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1235 IF (ibsign == 2) THEN
1236 wcart(1:3) = -wcart(1:3)
1237 kr = -kr
1238 END IF
1239 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1240 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1241 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1242
1243 found = .false.
1244 DO j = 1, nfull
1245 diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1246 diff(1:3) = diff(1:3) - anint(diff(1:3))
1247 IF (all(abs(diff(1:3)) < eps)) THEN
1248 found = .true.
1249 IF (.NOT. used(j)) THEN
1250 used(j) = .true.
1251 csym%kplink(1, j) = i
1252 csym%kpop(j) = kr
1253 wred(nred) = wred(nred) + 1.0_dp
1254 ELSE
1255 cpassert(csym%kplink(1, j) == i)
1256 END IF
1257 EXIT
1258 END IF
1259 END DO
1260 IF (.NOT. found) THEN
1261 CALL cp_abort(__location__, &
1262 "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1263 "to be closed under the K290 symmetry operations.")
1264 END IF
1265 END DO
1266 END DO
1267 END DO
1268
1269 CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1270
1271 DEALLOCATE (rep, used, wred)
1272
1273 END SUBROUTINE reduce_general_k290
1274
1275! **************************************************************************************************
1276!> \brief Reduce an explicit k-point set with SPGLIB symmetry operations.
1277!> \param csym ...
1278!> \param xkp_full explicit k-point coordinates
1279! **************************************************************************************************
1280 SUBROUTINE reduce_general_spglib(csym, xkp_full)
1281 TYPE(csym_type) :: csym
1282 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1283
1284 INTEGER :: i, iop, isign, j, kr, nfull, nred, nrot
1285 INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1286 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
1287 INTEGER, DIMENSION(3, 3) :: krot
1288 LOGICAL :: found
1289 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1290 REAL(kind=dp) :: eps
1291 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wred
1292 REAL(kind=dp), DIMENSION(3) :: diff, rr
1293
1294 nfull = SIZE(xkp_full, 2)
1295 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1296 ALLOCATE (srot(3, 3, csym%n_operations))
1297 CALL setup_spglib_operations(csym, srot, nrot)
1298
1299 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1300 used = .false.
1301 rep = 0
1302 wred = 0.0_dp
1303 nred = 0
1304
1305 DO i = 1, nfull
1306 IF (used(i)) cycle
1307 nred = nred + 1
1308 rep(nred) = i
1309 used(i) = .true.
1310 csym%kplink(1, i) = i
1311 csym%kpop(i) = 1
1312 wred(nred) = 1.0_dp
1313
1314 DO iop = 1, nrot
1315 kr = csym%ibrot(iop)
1316 krot = reciprocal_rotation(srot(:, :, kr))
1317 DO isign = 1, 2
1318 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp_full(1:3, i))
1319 IF (isign == 2) THEN
1320 rr(1:3) = -rr(1:3)
1321 kr = -csym%ibrot(iop)
1322 ELSE
1323 kr = csym%ibrot(iop)
1324 END IF
1325
1326 found = .false.
1327 DO j = 1, nfull
1328 diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1329 diff(1:3) = diff(1:3) - anint(diff(1:3))
1330 IF (all(abs(diff(1:3)) < eps)) THEN
1331 found = .true.
1332 IF (.NOT. used(j)) THEN
1333 used(j) = .true.
1334 csym%kplink(1, j) = i
1335 csym%kpop(j) = kr
1336 wred(nred) = wred(nred) + 1.0_dp
1337 ELSE
1338 cpassert(csym%kplink(1, j) == i)
1339 END IF
1340 EXIT
1341 END IF
1342 END DO
1343 IF (.NOT. found) THEN
1344 CALL cp_abort(__location__, &
1345 "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1346 "to be closed under the requested symmetry operations.")
1347 END IF
1348 END DO
1349 END DO
1350 END DO
1351
1352 CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1353
1354 DEALLOCATE (rep, srot, used, wred)
1355
1356 END SUBROUTINE reduce_general_spglib
1357
1358! **************************************************************************************************
1359!> \brief Reduce an explicit k-point set with SPGLIB rotations and K290 operations.
1360!> \param csym ...
1361!> \param xkp_full explicit k-point coordinates
1362! **************************************************************************************************
1363 SUBROUTINE reduce_general_spglib_k290(csym, xkp_full)
1364 TYPE(csym_type) :: csym
1365 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1366
1367 INTEGER :: i, iop, isign, j, k290_op, nfull, nred, &
1368 nrot, nskipped
1369 INTEGER, ALLOCATABLE, DIMENSION(:) :: rep
1370 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
1371 INTEGER, DIMENSION(3, 3) :: krot
1372 LOGICAL :: found, valid
1373 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
1374 REAL(kind=dp) :: alat, eps
1375 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wred
1376 REAL(kind=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, diff, rr
1377
1378 nfull = SIZE(xkp_full, 2)
1379 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1380 CALL setup_k290_lattice(csym, a1, a2, a3, b1, b2, b3, alat)
1381 ALLOCATE (srot(3, 3, csym%n_operations))
1382 CALL setup_spglib_reduction_rotations(csym, srot, nrot)
1383
1384 ALLOCATE (rep(nfull), used(nfull), wred(nfull))
1385 used = .false.
1386 rep = 0
1387 wred = 0.0_dp
1388 nred = 0
1389 nskipped = 0
1390
1391 DO i = 1, nfull
1392 IF (used(i)) cycle
1393 nred = nred + 1
1394 rep(nred) = i
1395 used(i) = .true.
1396 csym%kplink(1, i) = i
1397 csym%kpop(i) = 1
1398 wred(nred) = 1.0_dp
1399
1400 DO iop = 1, nrot
1401 krot = reciprocal_rotation(srot(:, :, iop))
1402 DO isign = 1, 2
1403 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp_full(1:3, i))
1404 IF (isign == 2) rr(1:3) = -rr(1:3)
1405
1406 found = .false.
1407 DO j = 1, nfull
1408 diff(1:3) = xkp_full(1:3, j) - rr(1:3)
1409 diff(1:3) = diff(1:3) - anint(diff(1:3))
1410 IF (all(abs(diff(1:3)) < eps)) THEN
1411 found = .true.
1412 CALL find_k290_kpoint_operation(csym, xkp_full(1:3, i), xkp_full(1:3, j), &
1413 a1, a2, a3, b1, b2, b3, alat, &
1414 k290_op, valid)
1415 IF (.NOT. valid) THEN
1416 nskipped = nskipped + 1
1417 EXIT
1418 END IF
1419 IF (.NOT. used(j)) THEN
1420 used(j) = .true.
1421 csym%kplink(1, j) = i
1422 csym%kpop(j) = k290_op
1423 wred(nred) = wred(nred) + 1.0_dp
1424 ELSE
1425 cpassert(csym%kplink(1, j) == i)
1426 END IF
1427 EXIT
1428 END IF
1429 END DO
1430 IF (.NOT. found) THEN
1431 CALL cp_abort(__location__, &
1432 "KPOINTS%SYMMETRY with SCHEME GENERAL requires the explicit k-point set "// &
1433 "to be closed under the SPGLIB symmetry operations.")
1434 END IF
1435 END DO
1436 END DO
1437 END DO
1438
1439 IF (nskipped > 0) THEN
1440 CALL cp_warn(__location__, &
1441 "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1442 "the GENERAL k-point set was reduced only by the compatible mappings.")
1443 END IF
1444
1445 CALL store_general_reduction(csym, xkp_full, rep, wred, nred)
1446
1447 DEALLOCATE (rep, srot, used, wred)
1448
1449 END SUBROUTINE reduce_general_spglib_k290
1450
1451! **************************************************************************************************
1452!> \brief Store reduced GENERAL k-point representatives.
1453!> \param csym ...
1454!> \param xkp_full explicit k-point coordinates
1455!> \param rep representative indices
1456!> \param wred representative multiplicities
1457!> \param nred number of reduced representatives
1458! **************************************************************************************************
1459 SUBROUTINE store_general_reduction(csym, xkp_full, rep, wred, nred)
1460 TYPE(csym_type) :: csym
1461 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp_full
1462 INTEGER, DIMENSION(:), INTENT(IN) :: rep
1463 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wred
1464 INTEGER, INTENT(IN) :: nred
1465
1466 INTEGER :: i, j, nfull
1467
1468 nfull = SIZE(xkp_full, 2)
1469 csym%nkpoint = nred
1470 ALLOCATE (csym%xkpoint(3, nred), csym%wkpoint(nred))
1471 DO i = 1, nred
1472 csym%xkpoint(1:3, i) = xkp_full(1:3, rep(i))
1473 csym%wkpoint(i) = wred(i)
1474 END DO
1475 DO i = 1, nfull
1476 DO j = 1, nred
1477 IF (csym%kplink(1, i) == rep(j)) THEN
1478 csym%kplink(2, i) = j
1479 EXIT
1480 END IF
1481 END DO
1482 cpassert(csym%kplink(2, i) /= 0)
1483 END DO
1484
1485 END SUBROUTINE store_general_reduction
1486
1487! **************************************************************************************************
1488!> \brief Reduce a k-point mesh with SPGLIB rotations and K290 operations
1489!> \param csym ...
1490!> \param xkp full k-point mesh in reciprocal lattice coordinates
1491!> \param wkp reduced k-point weights
1492!> \param kpop K290 operation mapping the representative k-point to a mesh point
1493!> \param srot SPGLIB integer rotations in fractional coordinates
1494!> \param nrot number of stored rotations
1495!> \param a1 first lattice vector
1496!> \param a2 second lattice vector
1497!> \param a3 third lattice vector
1498!> \param b1 first reciprocal lattice vector
1499!> \param b2 second reciprocal lattice vector
1500!> \param b3 third reciprocal lattice vector
1501!> \param alat lattice scaling used by K290
1502! **************************************************************************************************
1503 SUBROUTINE reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
1504 a1, a2, a3, b1, b2, b3, alat)
1505 TYPE(csym_type) :: csym
1506 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1507 REAL(kind=dp), DIMENSION(:) :: wkp
1508 INTEGER, DIMENSION(:) :: kpop
1509 INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
1510 INTEGER, INTENT(IN) :: nrot
1511 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1512 REAL(kind=dp), INTENT(IN) :: alat
1513
1514 INTEGER :: i, iop, isign, j, k290_op, nkpts, &
1515 nskipped
1516 INTEGER, DIMENSION(3, 3) :: krot
1517 LOGICAL :: valid
1518 REAL(kind=dp) :: eps
1519 REAL(kind=dp), DIMENSION(3) :: diff, rr
1520
1521 nkpts = SIZE(wkp)
1522 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1523 nskipped = 0
1524
1525 wkp = 0.0_dp
1526 kpop = 0
1527 csym%kplink(1, :) = 0
1528
1529 DO i = 1, nkpts
1530 IF (csym%kplink(1, i) /= 0) cycle
1531
1532 csym%kplink(1, i) = i
1533 wkp(i) = 1.0_dp
1534 kpop(i) = 1
1535
1536 DO iop = 1, nrot
1537 krot = reciprocal_rotation(srot(:, :, iop))
1538 DO isign = 1, 2
1539 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp(1:3, i))
1540 IF (isign == 2) rr(1:3) = -rr(1:3)
1541
1542 DO j = 1, nkpts
1543 diff(1:3) = xkp(1:3, j) - rr(1:3)
1544 diff(1:3) = diff(1:3) - anint(diff(1:3))
1545 IF (all(abs(diff(1:3)) < eps)) THEN
1546 IF (j == i) EXIT
1547 IF (csym%kplink(1, j) /= 0) THEN
1548 cpassert(csym%kplink(1, j) == i)
1549 EXIT
1550 END IF
1551
1552 CALL find_k290_kpoint_operation(csym, xkp(1:3, i), xkp(1:3, j), &
1553 a1, a2, a3, b1, b2, b3, alat, &
1554 k290_op, valid)
1555 IF (.NOT. valid) THEN
1556 nskipped = nskipped + 1
1557 EXIT
1558 END IF
1559 csym%kplink(1, j) = i
1560 wkp(i) = wkp(i) + 1.0_dp
1561 kpop(j) = k290_op
1562 EXIT
1563 END IF
1564 END DO
1565 IF (j > nkpts) cycle
1566 END DO
1567 END DO
1568 END DO
1569
1570 DO i = 1, nkpts
1571 cpassert(csym%kplink(1, i) /= 0)
1572 cpassert(kpop(i) /= 0)
1573 END DO
1574 IF (nskipped > 0) THEN
1575 CALL cp_warn(__location__, &
1576 "Some SPGLIB k-point mappings are not represented by the K290 backend; "// &
1577 "the mesh was reduced only by the compatible mappings.")
1578 END IF
1579
1580 END SUBROUTINE reduce_spglib_kpoint_mesh_k290
1581
1582! **************************************************************************************************
1583!> \brief Find a K290 operation that maps one fractional k-point to another
1584!> \param csym ...
1585!> \param xref representative k-point
1586!> \param xtarget target k-point
1587!> \param a1 first lattice vector
1588!> \param a2 second lattice vector
1589!> \param a3 third lattice vector
1590!> \param b1 first reciprocal lattice vector
1591!> \param b2 second reciprocal lattice vector
1592!> \param b3 third reciprocal lattice vector
1593!> \param alat lattice scaling used by K290
1594!> \param k290_op K290 operation identifier
1595!> \param valid whether a matching K290 operation was found
1596! **************************************************************************************************
1597 SUBROUTINE find_k290_kpoint_operation(csym, xref, xtarget, a1, a2, a3, b1, b2, b3, alat, &
1598 k290_op, valid)
1599 TYPE(csym_type) :: csym
1600 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xref, xtarget, a1, a2, a3, b1, b2, b3
1601 REAL(kind=dp), INTENT(IN) :: alat
1602 INTEGER, INTENT(OUT) :: k290_op
1603 LOGICAL, INTENT(OUT) :: valid
1604
1605 INTEGER :: ibsign, iop, kr
1606 REAL(kind=dp) :: eps
1607 REAL(kind=dp), DIMENSION(3) :: diff, rr, wcart
1608
1609 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1610 k290_op = 0
1611 valid = .false.
1612
1613 DO iop = 1, csym%nrtot
1614 IF (iop > SIZE(csym%rt, 3)) cycle
1615 IF (csym%ibrot(iop) == 0) cycle
1616 DO ibsign = 1, 2
1617 wcart(1:3) = alat*(xref(1)*b1(1:3) + xref(2)*b2(1:3) + xref(3)*b3(1:3))
1618 wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1619 IF (ibsign == 2) THEN
1620 wcart(1:3) = -wcart(1:3)
1621 kr = -csym%ibrot(iop)
1622 ELSE
1623 kr = csym%ibrot(iop)
1624 END IF
1625 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1626 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1627 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1628
1629 diff(1:3) = xtarget(1:3) - rr(1:3)
1630 diff(1:3) = diff(1:3) - anint(diff(1:3))
1631 IF (all(abs(diff(1:3)) < eps)) THEN
1632 k290_op = kr
1633 valid = .true.
1634 RETURN
1635 END IF
1636 END DO
1637 END DO
1638
1639 END SUBROUTINE find_k290_kpoint_operation
1640
1641! **************************************************************************************************
1642!> \brief Score SPGLIB operations to choose stable atom transformations
1643!> \param csym ...
1644!> \param iop operation index
1645!> \param srot integer rotation in fractional coordinates
1646!> \return score, lower values are preferred
1647! **************************************************************************************************
1648 FUNCTION spglib_operation_score(csym, iop, srot) RESULT(score)
1649 TYPE(csym_type), INTENT(IN) :: csym
1650 INTEGER, INTENT(IN) :: iop
1651 INTEGER, DIMENSION(3, 3), INTENT(IN) :: srot
1652 INTEGER :: score
1653
1654 INTEGER :: i, nat
1655 INTEGER, DIMENSION(3, 3) :: eye, r2
1656 REAL(kind=dp) :: eps
1657
1658 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1659 nat = SIZE(csym%f0, 1)
1660 score = 0
1661 DO i = 1, nat
1662 IF (csym%f0(i, iop) /= i) score = score + 100
1663 END DO
1664 IF (any(abs(csym%vt(1:3, iop) - anint(csym%vt(1:3, iop))) > eps)) score = score + 10
1665
1666 eye = 0
1667 eye(1, 1) = 1
1668 eye(2, 2) = 1
1669 eye(3, 3) = 1
1670 r2(1:3, 1:3) = matmul(srot(1:3, 1:3), srot(1:3, 1:3))
1671 IF (any(r2(1:3, 1:3) /= eye(1:3, 1:3))) score = score + 1
1672
1673 END FUNCTION spglib_operation_score
1674
1675! **************************************************************************************************
1676!> \brief Reciprocal-space rotation corresponding to a fractional direct-space rotation
1677!> \param rot direct-space rotation
1678!> \return reciprocal-space rotation
1679! **************************************************************************************************
1680 FUNCTION reciprocal_rotation(rot) RESULT(krot)
1681 INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
1682 INTEGER, DIMENSION(3, 3) :: krot
1683
1684 REAL(kind=dp), DIMENSION(3, 3) :: rinv
1685
1686 rinv = inv_3x3(real(rot(1:3, 1:3), kind=dp))
1687 krot(1:3, 1:3) = nint(transpose(rinv(1:3, 1:3)))
1688
1689 END FUNCTION reciprocal_rotation
1690
1691! **************************************************************************************************
1692!> \brief Reduce a CP2K Monkhorst-Pack mesh using K290 symmetry operations
1693!> \param csym ...
1694!> \param xkp full k-point mesh in reciprocal lattice coordinates
1695!> \param wkp reduced k-point weights
1696!> \param kpop symmetry operation mapping the representative k-point to a mesh point
1697!> \param nc number of point group operations
1698!> \param ib K290 operation identifiers
1699!> \param r K290 rotation matrices
1700!> \param a1 first lattice vector
1701!> \param a2 second lattice vector
1702!> \param a3 third lattice vector
1703!> \param b1 first reciprocal lattice vector
1704!> \param b2 second reciprocal lattice vector
1705!> \param b3 third reciprocal lattice vector
1706!> \param alat lattice scaling used by K290
1707! **************************************************************************************************
1708 SUBROUTINE reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
1709 TYPE(csym_type) :: csym
1710 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1711 REAL(kind=dp), DIMENSION(:) :: wkp
1712 INTEGER, DIMENSION(:) :: kpop
1713 INTEGER, INTENT(IN) :: nc
1714 INTEGER, DIMENSION(48), INTENT(IN) :: ib
1715 REAL(kind=dp), DIMENSION(3, 3, 48), INTENT(IN) :: r
1716 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1717 REAL(kind=dp), INTENT(IN) :: alat
1718
1719 INTEGER :: i, ibsign, iop, j, kr, nkpts
1720 REAL(kind=dp) :: eps
1721 REAL(kind=dp), DIMENSION(3) :: diff, rr, wcart
1722
1723 nkpts = SIZE(wkp)
1724 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
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 ELSE
1759 cpassert(csym%kplink(1, j) == i)
1760 END IF
1761 EXIT
1762 END IF
1763 END DO
1764 ! Some point-group operations are incompatible with the requested Monkhorst-Pack mesh.
1765 IF (j > nkpts) cycle
1766 END DO
1767 END DO
1768 END DO
1769
1770 DO i = 1, nkpts
1771 cpassert(csym%kplink(1, i) /= 0)
1772 cpassert(kpop(i) /= 0)
1773 END DO
1774
1775 END SUBROUTINE reduce_kpoint_mesh
1776
1777!> \brief ...
1778!> \param nk ...
1779!> \param xkp ...
1780!> \param wkp ...
1781!> \param shift ...
1782!> \param gamma_centered ...
1783! **************************************************************************************************
1784 SUBROUTINE full_grid_gen(nk, xkp, wkp, shift, gamma_centered)
1785 INTEGER, INTENT(IN) :: nk(3)
1786 REAL(kind=dp), DIMENSION(:, :) :: xkp
1787 REAL(kind=dp), DIMENSION(:) :: wkp
1788 REAL(kind=dp), INTENT(IN) :: shift(3)
1789 LOGICAL, INTENT(IN), OPTIONAL :: gamma_centered
1790
1791 INTEGER :: i, idim, ix, iy, iz
1792 INTEGER, DIMENSION(3) :: ik
1793 LOGICAL :: gamma_mesh
1794 REAL(kind=dp) :: kpt_latt(3)
1795
1796 IF (PRESENT(gamma_centered)) THEN
1797 gamma_mesh = gamma_centered
1798 ELSE
1799 gamma_mesh = .false.
1800 END IF
1801
1802 wkp = 0.0_dp
1803 i = 0
1804 DO ix = 1, nk(1)
1805 DO iy = 1, nk(2)
1806 DO iz = 1, nk(3)
1807 i = i + 1
1808 ik(1) = ix
1809 ik(2) = iy
1810 ik(3) = iz
1811 DO idim = 1, 3
1812 IF (gamma_mesh .AND. mod(nk(idim), 2) == 0) THEN
1813 kpt_latt(idim) = real(2*ik(idim) - nk(idim), kind=dp)/ &
1814 (2._dp*real(nk(idim), kind=dp))
1815 ELSE
1816 kpt_latt(idim) = real(2*ik(idim) - nk(idim) - 1, kind=dp)/ &
1817 (2._dp*real(nk(idim), kind=dp))
1818 END IF
1819 END DO
1820 xkp(1:3, i) = kpt_latt(1:3)
1821 wkp(i) = 1.0_dp
1822 END DO
1823 END DO
1824 END DO
1825 DO i = 1, nk(1)*nk(2)*nk(3)
1826 xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
1827 END DO
1828
1829 END SUBROUTINE full_grid_gen
1830
1831! **************************************************************************************************
1832!> \brief ...
1833!> \param xkp ...
1834!> \param wkp ...
1835!> \param link ...
1836! **************************************************************************************************
1837 SUBROUTINE inversion_symm(xkp, wkp, link)
1838 REAL(kind=dp), DIMENSION(:, :) :: xkp
1839 REAL(kind=dp), DIMENSION(:) :: wkp
1840 INTEGER, DIMENSION(:) :: link
1841
1842 INTEGER :: i, j, nkpts
1843 REAL(kind=dp), DIMENSION(3) :: diff
1844
1845 nkpts = SIZE(wkp, 1)
1846
1847 link(:) = 0
1848 DO i = 1, nkpts
1849 IF (link(i) == 0) link(i) = i
1850 DO j = i + 1, nkpts
1851 IF (wkp(j) == 0) cycle
1852 diff(1:3) = xkp(1:3, i) + xkp(1:3, j)
1853 diff(1:3) = diff(1:3) - anint(diff(1:3))
1854 IF (all(abs(diff(1:3)) < 1.e-12_dp)) THEN
1855 wkp(i) = wkp(i) + wkp(j)
1856 wkp(j) = 0.0_dp
1857 link(j) = i
1858 EXIT
1859 END IF
1860 END DO
1861 END DO
1862
1863 END SUBROUTINE inversion_symm
1864
1865! **************************************************************************************************
1866!> \brief ...
1867!> \param x ...
1868!> \param r ...
1869!> \return ...
1870! **************************************************************************************************
1871 FUNCTION kp_apply_operation(x, r) RESULT(y)
1872 REAL(kind=dp), INTENT(IN) :: x(3), r(3, 3)
1873 REAL(kind=dp) :: y(3)
1874
1875 y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
1876 y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
1877 y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
1878
1879 END FUNCTION kp_apply_operation
1880
1881! **************************************************************************************************
1882!> \brief ...
1883!> \param csym ...
1884! **************************************************************************************************
1885 SUBROUTINE print_crys_symmetry(csym)
1886 TYPE(csym_type) :: csym
1887
1888 INTEGER :: i, iunit, j, plevel
1889
1890 iunit = csym%punit
1891 IF (iunit >= 0) THEN
1892 plevel = csym%plevel
1893 WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
1894 IF (csym%symlib) THEN
1895 WRITE (iunit, '(A,T71,A10)') " International Symbol: ", adjustr(trim(csym%international_symbol))
1896 WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", adjustr(trim(csym%pointgroup_symbol))
1897 WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", adjustr(trim(csym%schoenflies))
1898 !
1899 WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
1900 IF (plevel > 0) THEN
1901 DO i = 1, csym%n_operations
1902 WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
1903 " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
1904 WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
1905 END DO
1906 END IF
1907 ELSE
1908 IF (csym%spglib_requested) THEN
1909 WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not available"
1910 ELSE
1911 WRITE (iunit, "(T2,A)") "SPGLIB Crystal Symmetry Information was not requested"
1912 END IF
1913 END IF
1914 END IF
1915
1916 END SUBROUTINE print_crys_symmetry
1917
1918! **************************************************************************************************
1919!> \brief ...
1920!> \param csym ...
1921! **************************************************************************************************
1922 SUBROUTINE print_kp_symmetry(csym)
1923 TYPE(csym_type), INTENT(IN) :: csym
1924
1925 INTEGER :: i, iunit, nat, nmesh, plevel
1926
1927 iunit = csym%punit
1928 IF (iunit >= 0) THEN
1929 plevel = csym%plevel
1930 WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
1931 WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
1932 WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
1933 DO i = 1, csym%nkpoint
1934 WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), nint(csym%wkpoint(i))
1935 END DO
1936 nmesh = csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
1937 IF (nmesh > 0) THEN
1938 WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
1939 ELSE
1940 nmesh = SIZE(csym%kpmesh, 2)
1941 WRITE (iunit, '(/,A,T70,I10)') " Explicit K-point Set: ", nmesh
1942 END IF
1943 WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
1944 DO i = 1, nmesh
1945 WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
1946 csym%kplink(1:2, i), csym%kpop(i)
1947 END DO
1948 IF (csym%nrtot > 0) THEN
1949 WRITE (iunit, '(/,A)') " Atom Transformation Table"
1950 nat = SIZE(csym%f0, 1)
1951 DO i = 1, csym%nrtot
1952 WRITE (iunit, '(T10,A,I5,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
1953 END DO
1954 END IF
1955 END IF
1956
1957 END SUBROUTINE print_kp_symmetry
1958
1959END 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:1886
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:1923
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