(git:495eafe)
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
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 INTEGER :: plevel = 0
51 INTEGER :: punit = -1
52 INTEGER :: istriz = -1
53 REAL(kind=dp) :: delta = 1.0e-8_dp
54 REAL(kind=dp), DIMENSION(3, 3) :: hmat = 0.0_dp
55 ! KPOINTS
56 REAL(kind=dp), DIMENSION(3) :: wvk0 = 0.0_dp
57 INTEGER, DIMENSION(3) :: mesh = 0
58 INTEGER :: nkpoint = 0
59 INTEGER :: nat = 0
60 INTEGER, DIMENSION(:), ALLOCATABLE :: atype
61 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: scoord
62 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: xkpoint
63 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: wkpoint
64 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: kpmesh
65 INTEGER, DIMENSION(:, :), ALLOCATABLE :: kplink
66 INTEGER, DIMENSION(:), ALLOCATABLE :: kpop
67 !SPGLIB
68 CHARACTER(len=11) :: international_symbol = ""
69 CHARACTER(len=6) :: pointgroup_symbol = ""
70 CHARACTER(len=10) :: schoenflies = ""
71 INTEGER :: n_operations = 0
72 INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: rotations
73 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: translations
74 !K290
75 REAL(kind=dp), DIMENSION(:, :, :), ALLOCATABLE :: rt
76 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: vt
77 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0
78 INTEGER :: nrtot = 0
79 INTEGER, DIMENSION(:), ALLOCATABLE :: ibrot
80 END TYPE csym_type
81
82CONTAINS
83
84! **************************************************************************************************
85!> \brief Release the CSYM type
86!> \param csym The CSYM type
87! **************************************************************************************************
88 SUBROUTINE release_csym_type(csym)
89 TYPE(csym_type) :: csym
90
91 IF (ALLOCATED(csym%rotations)) THEN
92 DEALLOCATE (csym%rotations)
93 END IF
94 IF (ALLOCATED(csym%translations)) THEN
95 DEALLOCATE (csym%translations)
96 END IF
97 IF (ALLOCATED(csym%atype)) THEN
98 DEALLOCATE (csym%atype)
99 END IF
100 IF (ALLOCATED(csym%scoord)) THEN
101 DEALLOCATE (csym%scoord)
102 END IF
103 IF (ALLOCATED(csym%xkpoint)) THEN
104 DEALLOCATE (csym%xkpoint)
105 END IF
106 IF (ALLOCATED(csym%wkpoint)) THEN
107 DEALLOCATE (csym%wkpoint)
108 END IF
109 IF (ALLOCATED(csym%kpmesh)) THEN
110 DEALLOCATE (csym%kpmesh)
111 END IF
112 IF (ALLOCATED(csym%kplink)) THEN
113 DEALLOCATE (csym%kplink)
114 END IF
115 IF (ALLOCATED(csym%kpop)) THEN
116 DEALLOCATE (csym%kpop)
117 END IF
118 IF (ALLOCATED(csym%rt)) THEN
119 DEALLOCATE (csym%rt)
120 END IF
121 IF (ALLOCATED(csym%vt)) THEN
122 DEALLOCATE (csym%vt)
123 END IF
124 IF (ALLOCATED(csym%f0)) THEN
125 DEALLOCATE (csym%f0)
126 END IF
127 IF (ALLOCATED(csym%ibrot)) THEN
128 DEALLOCATE (csym%ibrot)
129 END IF
130
131 END SUBROUTINE release_csym_type
132
133! **************************************************************************************************
134!> \brief ...
135!> \param csym ...
136!> \param scoor ...
137!> \param types ...
138!> \param hmat ...
139!> \param delta ...
140!> \param iounit ...
141! **************************************************************************************************
142 SUBROUTINE crys_sym_gen(csym, scoor, types, hmat, delta, iounit)
143 TYPE(csym_type) :: csym
144 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: scoor
145 INTEGER, DIMENSION(:), INTENT(IN) :: types
146 REAL(kind=dp), INTENT(IN) :: hmat(3, 3)
147 REAL(kind=dp), INTENT(IN), OPTIONAL :: delta
148 INTEGER, INTENT(IN), OPTIONAL :: iounit
149
150 CHARACTER(LEN=*), PARAMETER :: routinen = 'crys_sym_gen'
151
152 INTEGER :: handle, ierr, major, micro, minor, nat, &
153 nop, tra_mat(3, 3)
154 LOGICAL :: spglib
155
156 CALL timeset(routinen, handle)
157
158 !..total number of atoms
159 nat = SIZE(scoor, 2)
160 csym%nat = nat
161
162 ! output unit
163 IF (PRESENT(iounit)) THEN
164 csym%punit = iounit
165 ELSE
166 csym%punit = -1
167 END IF
168
169 ! accuracy for symmetry
170 IF (PRESENT(delta)) THEN
171 csym%delta = delta
172 ELSE
173 csym%delta = 1.e-6_dp
174 END IF
175
176 !..set cell values
177 csym%hmat = hmat
178
179 ! atom types
180 ALLOCATE (csym%atype(nat))
181 csym%atype(1:nat) = types(1:nat)
182
183 ! scaled coordinates
184 ALLOCATE (csym%scoord(3, nat))
185 csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
186
187 csym%n_operations = 0
188
189 !..try spglib
190 major = spg_get_major_version()
191 minor = spg_get_minor_version()
192 micro = spg_get_micro_version()
193 IF (major == 0) THEN
194 CALL cp_warn(__location__, "Symmetry library SPGLIB not available")
195 spglib = .false.
196 ELSE
197 spglib = .true.
198 CALL cite_reference(togo2018)
199 ierr = spg_get_international(csym%international_symbol, transpose(hmat), scoor, types, nat, delta)
200 IF (ierr == 0) THEN
201 CALL cp_warn(__location__, "Symmetry Library SPGLIB failed")
202 spglib = .false.
203 ELSE
204 nop = spg_get_multiplicity(transpose(hmat), scoor, types, nat, delta)
205 ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
206 csym%n_operations = nop
207 ierr = spg_get_symmetry(csym%rotations, csym%translations, nop, &
208 transpose(hmat), scoor, types, nat, delta)
209 ! Schoenflies Symbol
210 csym%schoenflies = ' '
211 ierr = spg_get_schoenflies(csym%schoenflies, transpose(hmat), scoor, types, nat, delta)
212 ! Point Group
213 csym%pointgroup_symbol = ' '
214 tra_mat = 0
215 ierr = spg_get_pointgroup(csym%pointgroup_symbol, tra_mat, &
216 csym%rotations, csym%n_operations)
217
218 CALL strip_control_codes(csym%international_symbol)
219 CALL strip_control_codes(csym%schoenflies)
220 CALL strip_control_codes(csym%pointgroup_symbol)
221 END IF
222 END IF
223 csym%symlib = spglib
224
225 CALL timestop(handle)
226
227 END SUBROUTINE crys_sym_gen
228
229! **************************************************************************************************
230!> \brief ...
231!> \param csym ...
232!> \param nk ...
233!> \param symm ...
234!> \param shift ...
235!> \param full_grid ...
236!> \param inversion_symmetry_only ...
237!> \param use_spglib_reduction ...
238!> \param use_spglib_backend ...
239! **************************************************************************************************
240 SUBROUTINE kpoint_gen(csym, nk, symm, shift, full_grid, inversion_symmetry_only, &
241 use_spglib_reduction, use_spglib_backend)
242 TYPE(csym_type) :: csym
243 INTEGER, INTENT(IN) :: nk(3)
244 LOGICAL, INTENT(IN), OPTIONAL :: symm
245 REAL(kind=dp), INTENT(IN), OPTIONAL :: shift(3)
246 LOGICAL, INTENT(IN), OPTIONAL :: full_grid, inversion_symmetry_only, &
247 use_spglib_reduction, &
248 use_spglib_backend
249
250 CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_gen'
251
252 INTEGER :: handle, i, ik, j, nkp, nkpts
253 INTEGER, ALLOCATABLE, DIMENSION(:) :: kpop, xptr
254 LOGICAL :: fullmesh, inversion_only, &
255 orthorhombic_cell, spglib_backend, &
256 spglib_reduction
257 REAL(kind=dp) :: cell_eps
258 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: wkp
259 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
260
261 CALL timeset(routinen, handle)
262
263 IF (PRESENT(shift)) THEN
264 csym%wvk0 = shift
265 ELSE
266 csym%wvk0 = 0.0_dp
267 END IF
268
269 csym%istriz = -1
270 IF (PRESENT(symm)) THEN
271 IF (symm) csym%istriz = 1
272 END IF
273
274 IF (PRESENT(full_grid)) THEN
275 fullmesh = full_grid
276 ELSE
277 fullmesh = .false.
278 END IF
279 csym%fullgrid = fullmesh
280
281 IF (PRESENT(inversion_symmetry_only)) THEN
282 inversion_only = inversion_symmetry_only
283 ELSE
284 inversion_only = .false.
285 END IF
286 csym%inversion_only = inversion_only
287
288 IF (PRESENT(use_spglib_reduction)) THEN
289 spglib_reduction = use_spglib_reduction
290 ELSE
291 spglib_reduction = .false.
292 END IF
293 csym%spglib_reduction = spglib_reduction
294
295 IF (PRESENT(use_spglib_backend)) THEN
296 spglib_backend = use_spglib_backend
297 ELSE
298 spglib_backend = .false.
299 END IF
300 csym%spglib_backend = spglib_backend
301
302 cell_eps = 10.0_dp*csym%delta*max(1.0_dp, maxval(abs(csym%hmat)))
303 orthorhombic_cell = abs(csym%hmat(1, 2)) + abs(csym%hmat(1, 3)) + &
304 abs(csym%hmat(2, 1)) + abs(csym%hmat(2, 3)) + &
305 abs(csym%hmat(3, 1)) + abs(csym%hmat(3, 2)) <= cell_eps
306
307 IF (spglib_backend .AND. .NOT. spglib_reduction) THEN
308 CALL cp_abort(__location__, &
309 "SYMMETRY_BACKEND SPGLIB requires SYMMETRY_REDUCTION_METHOD SPGLIB")
310 END IF
311
312 csym%nkpoint = 0
313 csym%mesh(1:3) = nk(1:3)
314 csym%nrtot = 0
315 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
316 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
317 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
318 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
319 ALLOCATE (csym%rt(3, 3, 0), csym%vt(3, 0), csym%ibrot(0), csym%f0(csym%nat, 0))
320
321 nkpts = nk(1)*nk(2)*nk(3)
322 ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
323 ! kp: link
324 ALLOCATE (csym%kplink(2, nkpts))
325 csym%kplink = 0
326 kpop = 0
327
328 ! go through all the options
329 IF (csym%symlib) THEN
330 ! symmetry library is available
331 IF (fullmesh) THEN
332 ! full mesh requested
333 CALL full_grid_gen(nk, xkp, wkp, shift)
334 IF (csym%istriz == 1) THEN
335 ! use inversion symmetry
336 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
337 ELSE
338 ! full kpoint mesh is used
339 END IF
340 ELSE IF (csym%istriz /= 1 .OR. inversion_only) THEN
341 ! use inversion symmetry
342 CALL full_grid_gen(nk, xkp, wkp, shift)
343 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
344 ELSE
345 ! use symmetry library to reduce k-points
346 IF (sum(abs(csym%wvk0)) /= 0.0_dp) THEN
347 CALL cp_abort(__location__, "MacDonald shifted k-point meshes are only "// &
348 "possible without symmetrization.")
349 END IF
350
351 CALL full_grid_gen(nk, xkp, wkp, shift)
352 IF (spglib_backend) THEN
353 CALL kp_symmetry_spglib(csym, xkp, wkp, kpop)
354 ELSE
355 CALL kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction=spglib_reduction)
356 END IF
357
358 END IF
359 ELSE
360 ! no symmetry library is available
361 CALL full_grid_gen(nk, xkp, wkp, shift)
362 IF (csym%istriz == 1 .AND. .NOT. fullmesh .AND. .NOT. inversion_only .AND. &
363 orthorhombic_cell) THEN
364 ! fall back to the K290 atom mapping for simple cells when SPGLIB is not linked
365 IF (sum(abs(csym%wvk0)) /= 0.0_dp) THEN
366 CALL cp_abort(__location__, "MacDonald shifted k-point meshes are only "// &
367 "possible without symmetrization.")
368 END IF
369 CALL kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction=.false.)
370 ELSE IF (csym%istriz /= 1 .AND. fullmesh) THEN
371 ! full kpoint mesh is used
372 DO i = 1, nkpts
373 csym%kplink(1, i) = i
374 END DO
375 ELSE
376 ! use inversion symmetry
377 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
378 END IF
379 END IF
380 ! count kpoints
381 nkp = 0
382 DO i = 1, nkpts
383 IF (wkp(i) > 0.0_dp) nkp = nkp + 1
384 END DO
385
386 ! store reduced kpoint set
387 csym%nkpoint = nkp
388 ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
389 ALLOCATE (xptr(nkp))
390 j = 0
391 DO ik = 1, nkpts
392 IF (wkp(ik) > 0.0_dp) THEN
393 j = j + 1
394 csym%wkpoint(j) = wkp(ik)
395 csym%xkpoint(1:3, j) = xkp(1:3, ik)
396 xptr(j) = ik
397 END IF
398 END DO
399 cpassert(j == nkp)
400
401 ! kp: mesh
402 ALLOCATE (csym%kpmesh(3, nkpts))
403 csym%kpmesh(1:3, 1:nkpts) = xkp(1:3, 1:nkpts)
404
405 ! kp: link
406 DO ik = 1, nkpts
407 i = csym%kplink(1, ik)
408 DO j = 1, nkp
409 IF (i == xptr(j)) THEN
410 csym%kplink(2, ik) = j
411 EXIT
412 END IF
413 END DO
414 END DO
415 DEALLOCATE (xptr)
416
417 ! kp: operations
418 ALLOCATE (csym%kpop(nkpts))
419 IF (csym%nrtot > 0 .AND. csym%istriz == 1 .AND. .NOT. fullmesh .AND. &
420 .NOT. inversion_only) THEN
421 ! atomic symmetry operations possible
422 csym%kpop(1:nkpts) = kpop(1:nkpts)
423 DO ik = 1, nkpts
424 cpassert(csym%kpop(ik) /= 0)
425 END DO
426 ELSE
427 ! only time reversal symmetry
428 DO ik = 1, nkpts
429 IF (wkp(ik) > 0.0_dp) THEN
430 csym%kpop(ik) = 1
431 ELSE
432 csym%kpop(ik) = 2
433 END IF
434 END DO
435 END IF
436
437 DEALLOCATE (xkp, wkp, kpop)
438
439 CALL timestop(handle)
440
441 END SUBROUTINE kpoint_gen
442
443! **************************************************************************************************
444!> \brief ...
445!> \param csym ...
446!> \param xkp ...
447!> \param wkp ...
448!> \param kpop ...
449!> \param use_spglib_reduction ...
450! **************************************************************************************************
451 SUBROUTINE kp_symmetry(csym, xkp, wkp, kpop, use_spglib_reduction)
452 TYPE(csym_type) :: csym
453 REAL(kind=dp), DIMENSION(:, :) :: xkp
454 REAL(kind=dp), DIMENSION(:) :: wkp
455 INTEGER, DIMENSION(:) :: kpop
456 LOGICAL, INTENT(IN), OPTIONAL :: use_spglib_reduction
457
458 INTEGER :: i, ihc, ihg, indpg, iou, iq1, iq2, iq3, &
459 istriz, isy, li, nat, nc, nhash, &
460 nkpoint, nrot, nsp, ntvec
461 INTEGER, ALLOCATABLE, DIMENSION(:) :: includ, isc, list, lwght, ty
462 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0, lrot
463 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
464 INTEGER, DIMENSION(48) :: ib
465 LOGICAL :: spglib_reduction
466 REAL(kind=dp) :: alat
467 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rlist, rx, tvec, wvkl, xkapa
468 REAL(kind=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, origin, wvk0
469 REAL(kind=dp), DIMENSION(3, 3) :: hmat, strain
470 REAL(kind=dp), DIMENSION(3, 3, 48) :: r
471 REAL(kind=dp), DIMENSION(3, 48) :: vt
472
473 iou = csym%punit
474 hmat = csym%hmat
475 nat = csym%nat
476 iq1 = csym%mesh(1)
477 iq2 = csym%mesh(2)
478 iq3 = csym%mesh(3)
479 nkpoint = 10*iq1*iq2*iq3
480 wvk0 = csym%wvk0
481 istriz = csym%istriz
482 IF (PRESENT(use_spglib_reduction)) THEN
483 spglib_reduction = use_spglib_reduction
484 ELSE
485 spglib_reduction = .false.
486 END IF
487 a1(1:3) = hmat(1:3, 1)
488 a2(1:3) = hmat(1:3, 2)
489 a3(1:3) = hmat(1:3, 3)
490 alat = hmat(1, 1)
491 strain = 0.0_dp
492 ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, 200), ty(nat), isc(nat), f0(49, nat))
493 ty(1:nat) = csym%atype(1:nat)
494 nsp = maxval(ty)
495 DO i = 1, nat
496 xkapa(1:3, i) = matmul(hmat, csym%scoord(1:3, i))
497 END DO
498 nhash = 1000
499 ALLOCATE (wvkl(3, nkpoint), rlist(3, nkpoint), includ(nkpoint), list(nhash + nkpoint))
500 ALLOCATE (lrot(48, nkpoint), lwght(nkpoint))
501
502 IF (iou > 0) THEN
503 WRITE (iou, '(/,(T2,A79))') &
504 "*******************************************************************************", &
505 "** Special K-Point Generation by K290 **", &
506 "*******************************************************************************"
507 END IF
508 CALL cite_reference(worlton1972)
509 IF (spglib_reduction) CALL cite_reference(togo2018)
510
511 CALL k290s(iou, nat, nkpoint, nsp, iq1, iq2, iq3, istriz, &
512 a1, a2, a3, alat, strain, xkapa, rx, tvec, &
513 ty, isc, f0, ntvec, wvk0, wvkl, lwght, lrot, &
514 nhash, includ, list, rlist, csym%delta)
515
516 CALL group1s(0, a1, a2, a3, nat, ty, xkapa, b1, b2, b3, &
517 ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
518 vt, f0, r, tvec, origin, rx, isc, csym%delta)
519
520 IF (iou > 0) THEN
521 WRITE (iou, '((T2,A79))') &
522 "*******************************************************************************", &
523 "** Finished K290 **", &
524 "*******************************************************************************"
525 END IF
526
527 csym%nrtot = nc
528 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
529 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
530 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
531 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
532 ALLOCATE (csym%rt(3, 3, nc), csym%vt(3, nc), csym%ibrot(nc))
533 csym%vt(1:3, 1:nc) = vt(1:3, 1:nc)
534 ALLOCATE (csym%f0(nat, nc))
535 DO i = 1, nc
536 csym%rt(1:3, 1:3, i) = r(1:3, 1:3, ib(i))
537 csym%f0(1:nat, i) = f0(i, 1:nat)
538 END DO
539 csym%ibrot(1:nc) = ib(1:nc)
540
541 IF (csym%n_operations > nc) THEN
542 IF (ALLOCATED(srot)) DEALLOCATE (srot)
543 ALLOCATE (srot(3, 3, csym%n_operations))
544 CALL setup_spglib_operations(csym, srot, nrot)
545 CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
546 DEALLOCATE (srot)
547 ELSE IF (spglib_reduction) THEN
548 ALLOCATE (srot(3, 3, csym%n_operations))
549 CALL setup_spglib_reduction_rotations(csym, srot, nrot)
550 CALL reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
551 a1, a2, a3, b1, b2, b3, alat)
552 DEALLOCATE (srot)
553 ELSE
554 CALL reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
555 END IF
556 DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
557 DEALLOCATE (wvkl, rlist, includ, list)
558 DEALLOCATE (lrot, lwght)
559
560 END SUBROUTINE kp_symmetry
561
562! **************************************************************************************************
563!> \brief Reduce a CP2K Monkhorst-Pack mesh using SPGLIB symmetry operations
564!> \param csym ...
565!> \param xkp ...
566!> \param wkp ...
567!> \param kpop ...
568! **************************************************************************************************
569 SUBROUTINE kp_symmetry_spglib(csym, xkp, wkp, kpop)
570 TYPE(csym_type) :: csym
571 REAL(kind=dp), DIMENSION(:, :) :: xkp
572 REAL(kind=dp), DIMENSION(:) :: wkp
573 INTEGER, DIMENSION(:) :: kpop
574
575 INTEGER :: iou, nrot
576 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: srot
577
578 iou = csym%punit
579 IF (iou > 0) THEN
580 WRITE (iou, '(/,(T2,A79))') &
581 "*******************************************************************************", &
582 "** Special K-Point Generation by SPGLIB **", &
583 "*******************************************************************************"
584 END IF
585 CALL cite_reference(togo2018)
586
587 ALLOCATE (srot(3, 3, csym%n_operations))
588 CALL setup_spglib_operations(csym, srot, nrot)
589 CALL reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
590 DEALLOCATE (srot)
591
592 IF (iou > 0) THEN
593 WRITE (iou, '((T2,A79))') &
594 "*******************************************************************************", &
595 "** Finished SPGLIB **", &
596 "*******************************************************************************"
597 END IF
598
599 END SUBROUTINE kp_symmetry_spglib
600
601! **************************************************************************************************
602!> \brief Store usable SPGLIB space-group operations for k-point symmetry
603!> \param csym ...
604!> \param srot integer rotations in fractional coordinates
605!> \param nrot number of stored rotations
606! **************************************************************************************************
607 SUBROUTINE setup_spglib_operations(csym, srot, nrot)
608 TYPE(csym_type) :: csym
609 INTEGER, DIMENSION(:, :, :), INTENT(OUT) :: srot
610 INTEGER, INTENT(OUT) :: nrot
611
612 INTEGER :: iop, jop, pass
613 INTEGER, ALLOCATABLE, DIMENSION(:) :: perm
614 INTEGER, DIMENSION(3, 3) :: eye, irot
615 LOGICAL :: duplicate, identity, valid, &
616 zero_translation
617 REAL(kind=dp) :: eps
618 REAL(kind=dp), DIMENSION(3, 3) :: h_inv, rfrac
619
620 cpassert(csym%symlib)
621
622 srot = 0
623 csym%nrtot = 0
624 IF (ALLOCATED(csym%rt)) DEALLOCATE (csym%rt)
625 IF (ALLOCATED(csym%vt)) DEALLOCATE (csym%vt)
626 IF (ALLOCATED(csym%ibrot)) DEALLOCATE (csym%ibrot)
627 IF (ALLOCATED(csym%f0)) DEALLOCATE (csym%f0)
628 ALLOCATE (csym%rt(3, 3, csym%n_operations), csym%vt(3, csym%n_operations))
629 ALLOCATE (csym%ibrot(csym%n_operations), csym%f0(csym%nat, csym%n_operations))
630 csym%rt = 0.0_dp
631 csym%vt = 0.0_dp
632 csym%ibrot = 0
633 csym%f0 = 0
634
635 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
636 h_inv = inv_3x3(csym%hmat)
637 ALLOCATE (perm(csym%nat))
638
639 eye = 0
640 eye(1, 1) = 1
641 eye(2, 2) = 1
642 eye(3, 3) = 1
643
644 nrot = 0
645 ! Operation 1 is used as the untransformed representative k-point.
646 ! Prefer integer translations before fractional alternatives with the same rotation.
647 DO pass = 1, 4
648 DO iop = 1, csym%n_operations
649 irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
650 identity = all(irot == eye)
651 zero_translation = all(abs(csym%translations(1:3, iop) - &
652 anint(csym%translations(1:3, iop))) <= eps)
653 IF (pass == 1 .AND. (.NOT. identity .OR. .NOT. zero_translation)) cycle
654 IF (pass == 2 .AND. (identity .OR. .NOT. zero_translation)) cycle
655 IF (pass == 3 .AND. (.NOT. identity .OR. zero_translation)) cycle
656 IF (pass == 4 .AND. (identity .OR. zero_translation)) cycle
657
658 duplicate = .false.
659 DO jop = 1, nrot
660 IF (all(irot == srot(:, :, jop)) .AND. &
661 all(abs(csym%translations(1:3, iop) - csym%vt(1:3, jop) - &
662 anint(csym%translations(1:3, iop) - csym%vt(1:3, jop))) <= eps)) THEN
663 duplicate = .true.
664 EXIT
665 END IF
666 END DO
667 IF (duplicate) cycle
668
669 CALL spglib_atom_permutation(csym, irot, csym%translations(:, iop), perm, valid)
670 IF (.NOT. valid) cycle
671
672 nrot = nrot + 1
673
674 srot(1:3, 1:3, nrot) = irot(1:3, 1:3)
675 rfrac(1:3, 1:3) = real(irot(1:3, 1:3), kind=dp)
676 csym%rt(1:3, 1:3, nrot) = matmul(csym%hmat, matmul(rfrac, h_inv))
677 csym%vt(1:3, nrot) = csym%translations(1:3, iop)
678 csym%ibrot(nrot) = nrot
679 csym%f0(1:csym%nat, nrot) = perm(1:csym%nat)
680 END DO
681 END DO
682
683 DEALLOCATE (perm)
684 csym%nrtot = nrot
685 IF (nrot == 0) CALL cp_abort(__location__, "SPGLIB did not return usable symmetry operations")
686
687 END SUBROUTINE setup_spglib_operations
688
689! **************************************************************************************************
690!> \brief Store unique SPGLIB rotations for K290-backend diagnostic reduction
691!> \param csym ...
692!> \param srot integer rotations in fractional coordinates
693!> \param nrot number of stored rotations
694! **************************************************************************************************
695 SUBROUTINE setup_spglib_reduction_rotations(csym, srot, nrot)
696 TYPE(csym_type) :: csym
697 INTEGER, DIMENSION(:, :, :), INTENT(OUT) :: srot
698 INTEGER, INTENT(OUT) :: nrot
699
700 INTEGER :: iop, jop, pass
701 INTEGER, DIMENSION(3, 3) :: eye, irot
702 LOGICAL :: duplicate, identity
703
704 cpassert(csym%symlib)
705
706 srot = 0
707 eye = 0
708 eye(1, 1) = 1
709 eye(2, 2) = 1
710 eye(3, 3) = 1
711
712 nrot = 0
713 ! Keep the identity first, matching the representative k-point operation.
714 DO pass = 1, 2
715 DO iop = 1, csym%n_operations
716 irot(1:3, 1:3) = csym%rotations(1:3, 1:3, iop)
717 identity = all(irot == eye)
718 IF (pass == 1 .AND. .NOT. identity) cycle
719 IF (pass == 2 .AND. identity) cycle
720
721 duplicate = .false.
722 DO jop = 1, nrot
723 IF (all(irot == srot(:, :, jop))) THEN
724 duplicate = .true.
725 EXIT
726 END IF
727 END DO
728 IF (duplicate) cycle
729
730 nrot = nrot + 1
731 srot(1:3, 1:3, nrot) = irot(1:3, 1:3)
732 END DO
733 END DO
734
735 IF (nrot == 0) CALL cp_abort(__location__, "SPGLIB did not return usable symmetry rotations")
736
737 END SUBROUTINE setup_spglib_reduction_rotations
738
739! **************************************************************************************************
740!> \brief Determine the atom permutation generated by a SPGLIB space-group operation
741!> \param csym ...
742!> \param rot integer rotation in fractional coordinates
743!> \param trans fractional translation
744!> \param perm atom permutation
745!> \param valid whether all atoms were mapped
746! **************************************************************************************************
747 SUBROUTINE spglib_atom_permutation(csym, rot, trans, perm, valid)
748 TYPE(csym_type) :: csym
749 INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
750 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: trans
751 INTEGER, DIMENSION(:), INTENT(OUT) :: perm
752 LOGICAL, INTENT(OUT) :: valid
753
754 INTEGER :: i, j, nat
755 LOGICAL :: found
756 LOGICAL, ALLOCATABLE, DIMENSION(:) :: used
757 REAL(kind=dp) :: eps
758 REAL(kind=dp), DIMENSION(3) :: diff, spos
759 REAL(kind=dp), DIMENSION(3, 3) :: rfrac
760
761 nat = csym%nat
762 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
763 rfrac(1:3, 1:3) = real(rot(1:3, 1:3), kind=dp)
764 ALLOCATE (used(nat))
765 used = .false.
766 perm = 0
767 valid = .true.
768
769 DO i = 1, nat
770 spos(1:3) = matmul(rfrac(1:3, 1:3), csym%scoord(1:3, i)) + trans(1:3)
771 found = .false.
772 DO j = 1, nat
773 IF (used(j)) cycle
774 IF (csym%atype(i) /= csym%atype(j)) cycle
775 diff(1:3) = spos(1:3) - csym%scoord(1:3, j)
776 diff(1:3) = diff(1:3) - anint(diff(1:3))
777 IF (all(abs(diff(1:3)) < eps)) THEN
778 perm(i) = j
779 used(j) = .true.
780 found = .true.
781 EXIT
782 END IF
783 END DO
784 IF (.NOT. found) THEN
785 valid = .false.
786 EXIT
787 END IF
788 END DO
789
790 DEALLOCATE (used)
791
792 END SUBROUTINE spglib_atom_permutation
793
794! **************************************************************************************************
795!> \brief Reduce a k-point mesh with SPGLIB direct-space operations
796!> \param csym ...
797!> \param xkp full k-point mesh in reciprocal lattice coordinates
798!> \param wkp reduced k-point weights
799!> \param kpop symmetry operation mapping the representative k-point to a mesh point
800!> \param srot integer rotations in fractional coordinates
801!> \param nrot number of stored rotations
802! **************************************************************************************************
803 SUBROUTINE reduce_spglib_kpoint_mesh(csym, xkp, wkp, kpop, srot, nrot)
804 TYPE(csym_type) :: csym
805 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
806 REAL(kind=dp), DIMENSION(:) :: wkp
807 INTEGER, DIMENSION(:) :: kpop
808 INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
809 INTEGER, INTENT(IN) :: nrot
810
811 INTEGER :: i, iop, isign, j, kr, nkpts, score
812 INTEGER, ALLOCATABLE, DIMENSION(:) :: kscore
813 INTEGER, DIMENSION(3, 3) :: krot
814 REAL(kind=dp) :: eps
815 REAL(kind=dp), DIMENSION(3) :: diff, rr
816
817 nkpts = SIZE(wkp)
818 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
819 ALLOCATE (kscore(nkpts))
820
821 wkp = 0.0_dp
822 kpop = 0
823 csym%kplink(1, :) = 0
824 kscore = huge(0)
825
826 DO i = 1, nkpts
827 IF (csym%kplink(1, i) /= 0) cycle
828
829 csym%kplink(1, i) = i
830 wkp(i) = 1.0_dp
831 kpop(i) = 1
832 kscore(i) = 0
833
834 DO iop = 1, nrot
835 kr = csym%ibrot(iop)
836 krot = reciprocal_rotation(srot(:, :, kr))
837 score = spglib_operation_score(csym, iop, srot(:, :, kr))
838 DO isign = 1, 2
839 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp(1:3, i))
840 IF (isign == 2) THEN
841 rr(1:3) = -rr(1:3)
842 kr = -csym%ibrot(iop)
843 ELSE
844 kr = csym%ibrot(iop)
845 END IF
846
847 DO j = 1, nkpts
848 diff(1:3) = xkp(1:3, j) - rr(1:3)
849 diff(1:3) = diff(1:3) - anint(diff(1:3))
850 IF (all(abs(diff(1:3)) < eps)) THEN
851 IF (csym%kplink(1, j) == 0) THEN
852 csym%kplink(1, j) = i
853 wkp(i) = wkp(i) + 1.0_dp
854 kpop(j) = kr
855 kscore(j) = score
856 ELSE
857 cpassert(csym%kplink(1, j) == i)
858 IF (score < kscore(j)) THEN
859 kpop(j) = kr
860 kscore(j) = score
861 END IF
862 END IF
863 EXIT
864 END IF
865 END DO
866 IF (j > nkpts) cycle
867 END DO
868 END DO
869 END DO
870
871 DO i = 1, nkpts
872 cpassert(csym%kplink(1, i) /= 0)
873 cpassert(kpop(i) /= 0)
874 END DO
875 DEALLOCATE (kscore)
876
877 END SUBROUTINE reduce_spglib_kpoint_mesh
878
879! **************************************************************************************************
880!> \brief Reduce a k-point mesh with SPGLIB rotations and K290 operations
881!> \param csym ...
882!> \param xkp full k-point mesh in reciprocal lattice coordinates
883!> \param wkp reduced k-point weights
884!> \param kpop K290 operation mapping the representative k-point to a mesh point
885!> \param srot SPGLIB integer rotations in fractional coordinates
886!> \param nrot number of stored rotations
887!> \param a1 first lattice vector
888!> \param a2 second lattice vector
889!> \param a3 third lattice vector
890!> \param b1 first reciprocal lattice vector
891!> \param b2 second reciprocal lattice vector
892!> \param b3 third reciprocal lattice vector
893!> \param alat lattice scaling used by K290
894! **************************************************************************************************
895 SUBROUTINE reduce_spglib_kpoint_mesh_k290(csym, xkp, wkp, kpop, srot, nrot, &
896 a1, a2, a3, b1, b2, b3, alat)
897 TYPE(csym_type) :: csym
898 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
899 REAL(kind=dp), DIMENSION(:) :: wkp
900 INTEGER, DIMENSION(:) :: kpop
901 INTEGER, DIMENSION(:, :, :), INTENT(IN) :: srot
902 INTEGER, INTENT(IN) :: nrot
903 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
904 REAL(kind=dp), INTENT(IN) :: alat
905
906 INTEGER :: i, iop, isign, j, k290_op, nkpts
907 INTEGER, DIMENSION(3, 3) :: krot
908 LOGICAL :: valid
909 REAL(kind=dp) :: eps
910 REAL(kind=dp), DIMENSION(3) :: diff, rr
911
912 nkpts = SIZE(wkp)
913 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
914
915 wkp = 0.0_dp
916 kpop = 0
917 csym%kplink(1, :) = 0
918
919 DO i = 1, nkpts
920 IF (csym%kplink(1, i) /= 0) cycle
921
922 csym%kplink(1, i) = i
923 wkp(i) = 1.0_dp
924 kpop(i) = 1
925
926 DO iop = 1, nrot
927 krot = reciprocal_rotation(srot(:, :, iop))
928 DO isign = 1, 2
929 rr(1:3) = matmul(real(krot(1:3, 1:3), kind=dp), xkp(1:3, i))
930 IF (isign == 2) rr(1:3) = -rr(1:3)
931
932 DO j = 1, nkpts
933 diff(1:3) = xkp(1:3, j) - rr(1:3)
934 diff(1:3) = diff(1:3) - anint(diff(1:3))
935 IF (all(abs(diff(1:3)) < eps)) THEN
936 IF (j == i) EXIT
937 IF (csym%kplink(1, j) /= 0) THEN
938 cpassert(csym%kplink(1, j) == i)
939 EXIT
940 END IF
941
942 CALL find_k290_kpoint_operation(csym, xkp(1:3, i), xkp(1:3, j), &
943 a1, a2, a3, b1, b2, b3, alat, &
944 k290_op, valid)
945 IF (.NOT. valid) THEN
946 CALL cp_abort(__location__, &
947 "SPGLIB k-point mapping is not represented by K290 backend")
948 END IF
949 csym%kplink(1, j) = i
950 wkp(i) = wkp(i) + 1.0_dp
951 kpop(j) = k290_op
952 EXIT
953 END IF
954 END DO
955 IF (j > nkpts) cycle
956 END DO
957 END DO
958 END DO
959
960 DO i = 1, nkpts
961 cpassert(csym%kplink(1, i) /= 0)
962 cpassert(kpop(i) /= 0)
963 END DO
964
965 END SUBROUTINE reduce_spglib_kpoint_mesh_k290
966
967! **************************************************************************************************
968!> \brief Find a K290 operation that maps one fractional k-point to another
969!> \param csym ...
970!> \param xref representative k-point
971!> \param xtarget target k-point
972!> \param a1 first lattice vector
973!> \param a2 second lattice vector
974!> \param a3 third lattice vector
975!> \param b1 first reciprocal lattice vector
976!> \param b2 second reciprocal lattice vector
977!> \param b3 third reciprocal lattice vector
978!> \param alat lattice scaling used by K290
979!> \param k290_op K290 operation identifier
980!> \param valid whether a matching K290 operation was found
981! **************************************************************************************************
982 SUBROUTINE find_k290_kpoint_operation(csym, xref, xtarget, a1, a2, a3, b1, b2, b3, alat, &
983 k290_op, valid)
984 TYPE(csym_type) :: csym
985 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: xref, xtarget, a1, a2, a3, b1, b2, b3
986 REAL(kind=dp), INTENT(IN) :: alat
987 INTEGER, INTENT(OUT) :: k290_op
988 LOGICAL, INTENT(OUT) :: valid
989
990 INTEGER :: ibsign, iop, kr
991 REAL(kind=dp) :: eps
992 REAL(kind=dp), DIMENSION(3) :: diff, rr, wcart
993
994 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
995 k290_op = 0
996 valid = .false.
997
998 DO iop = 1, csym%nrtot
999 IF (iop > SIZE(csym%rt, 3)) cycle
1000 IF (csym%ibrot(iop) == 0) cycle
1001 DO ibsign = 1, 2
1002 wcart(1:3) = alat*(xref(1)*b1(1:3) + xref(2)*b2(1:3) + xref(3)*b3(1:3))
1003 wcart(1:3) = kp_apply_operation(wcart(1:3), csym%rt(1:3, 1:3, iop))
1004 IF (ibsign == 2) THEN
1005 wcart(1:3) = -wcart(1:3)
1006 kr = -csym%ibrot(iop)
1007 ELSE
1008 kr = csym%ibrot(iop)
1009 END IF
1010 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1011 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1012 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1013
1014 diff(1:3) = xtarget(1:3) - rr(1:3)
1015 diff(1:3) = diff(1:3) - anint(diff(1:3))
1016 IF (all(abs(diff(1:3)) < eps)) THEN
1017 k290_op = kr
1018 valid = .true.
1019 RETURN
1020 END IF
1021 END DO
1022 END DO
1023
1024 END SUBROUTINE find_k290_kpoint_operation
1025
1026! **************************************************************************************************
1027!> \brief Score SPGLIB operations to choose stable atom transformations
1028!> \param csym ...
1029!> \param iop operation index
1030!> \param srot integer rotation in fractional coordinates
1031!> \return score, lower values are preferred
1032! **************************************************************************************************
1033 FUNCTION spglib_operation_score(csym, iop, srot) RESULT(score)
1034 TYPE(csym_type), INTENT(IN) :: csym
1035 INTEGER, INTENT(IN) :: iop
1036 INTEGER, DIMENSION(3, 3), INTENT(IN) :: srot
1037 INTEGER :: score
1038
1039 INTEGER :: i, nat
1040 INTEGER, DIMENSION(3, 3) :: eye, r2
1041 REAL(kind=dp) :: eps
1042
1043 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1044 nat = SIZE(csym%f0, 1)
1045 score = 0
1046 DO i = 1, nat
1047 IF (csym%f0(i, iop) /= i) score = score + 100
1048 END DO
1049 IF (any(abs(csym%vt(1:3, iop) - anint(csym%vt(1:3, iop))) > eps)) score = score + 10
1050
1051 eye = 0
1052 eye(1, 1) = 1
1053 eye(2, 2) = 1
1054 eye(3, 3) = 1
1055 r2(1:3, 1:3) = matmul(srot(1:3, 1:3), srot(1:3, 1:3))
1056 IF (any(r2(1:3, 1:3) /= eye(1:3, 1:3))) score = score + 1
1057
1058 END FUNCTION spglib_operation_score
1059
1060! **************************************************************************************************
1061!> \brief Reciprocal-space rotation corresponding to a fractional direct-space rotation
1062!> \param rot direct-space rotation
1063!> \return reciprocal-space rotation
1064! **************************************************************************************************
1065 FUNCTION reciprocal_rotation(rot) RESULT(krot)
1066 INTEGER, DIMENSION(3, 3), INTENT(IN) :: rot
1067 INTEGER, DIMENSION(3, 3) :: krot
1068
1069 REAL(kind=dp), DIMENSION(3, 3) :: rinv
1070
1071 rinv = inv_3x3(real(rot(1:3, 1:3), kind=dp))
1072 krot(1:3, 1:3) = nint(transpose(rinv(1:3, 1:3)))
1073
1074 END FUNCTION reciprocal_rotation
1075
1076! **************************************************************************************************
1077!> \brief Reduce a CP2K Monkhorst-Pack mesh using K290 symmetry operations
1078!> \param csym ...
1079!> \param xkp full k-point mesh in reciprocal lattice coordinates
1080!> \param wkp reduced k-point weights
1081!> \param kpop symmetry operation mapping the representative k-point to a mesh point
1082!> \param nc number of point group operations
1083!> \param ib K290 operation identifiers
1084!> \param r K290 rotation matrices
1085!> \param a1 first lattice vector
1086!> \param a2 second lattice vector
1087!> \param a3 third lattice vector
1088!> \param b1 first reciprocal lattice vector
1089!> \param b2 second reciprocal lattice vector
1090!> \param b3 third reciprocal lattice vector
1091!> \param alat lattice scaling used by K290
1092! **************************************************************************************************
1093 SUBROUTINE reduce_kpoint_mesh(csym, xkp, wkp, kpop, nc, ib, r, a1, a2, a3, b1, b2, b3, alat)
1094 TYPE(csym_type) :: csym
1095 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: xkp
1096 REAL(kind=dp), DIMENSION(:) :: wkp
1097 INTEGER, DIMENSION(:) :: kpop
1098 INTEGER, INTENT(IN) :: nc
1099 INTEGER, DIMENSION(48), INTENT(IN) :: ib
1100 REAL(kind=dp), DIMENSION(3, 3, 48), INTENT(IN) :: r
1101 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3, b1, b2, b3
1102 REAL(kind=dp), INTENT(IN) :: alat
1103
1104 INTEGER :: i, ibsign, iop, j, kr, nkpts
1105 REAL(kind=dp) :: eps
1106 REAL(kind=dp), DIMENSION(3) :: diff, rr, wcart
1107
1108 nkpts = SIZE(wkp)
1109 eps = max(1.e-12_dp, 10.0_dp*csym%delta)
1110
1111 wkp = 0.0_dp
1112 kpop = 0
1113 csym%kplink(1, :) = 0
1114
1115 DO i = 1, nkpts
1116 IF (csym%kplink(1, i) /= 0) cycle
1117
1118 csym%kplink(1, i) = i
1119 wkp(i) = 1.0_dp
1120 kpop(i) = 1
1121
1122 DO iop = 1, nc
1123 DO ibsign = 1, 2
1124 kr = ib(iop)
1125 wcart(1:3) = alat*(xkp(1, i)*b1(1:3) + xkp(2, i)*b2(1:3) + xkp(3, i)*b3(1:3))
1126 wcart(1:3) = kp_apply_operation(wcart(1:3), r(1:3, 1:3, kr))
1127 IF (ibsign == 2) THEN
1128 wcart(1:3) = -wcart(1:3)
1129 kr = -kr
1130 END IF
1131 rr(1) = dot_product(a1(1:3), wcart(1:3))/alat
1132 rr(2) = dot_product(a2(1:3), wcart(1:3))/alat
1133 rr(3) = dot_product(a3(1:3), wcart(1:3))/alat
1134
1135 DO j = 1, nkpts
1136 diff(1:3) = xkp(1:3, j) - rr(1:3)
1137 diff(1:3) = diff(1:3) - anint(diff(1:3))
1138 IF (all(abs(diff(1:3)) < eps)) THEN
1139 IF (csym%kplink(1, j) == 0) THEN
1140 csym%kplink(1, j) = i
1141 wkp(i) = wkp(i) + 1.0_dp
1142 kpop(j) = kr
1143 ELSE
1144 cpassert(csym%kplink(1, j) == i)
1145 END IF
1146 EXIT
1147 END IF
1148 END DO
1149 ! Some point-group operations are incompatible with the requested Monkhorst-Pack mesh.
1150 IF (j > nkpts) cycle
1151 END DO
1152 END DO
1153 END DO
1154
1155 DO i = 1, nkpts
1156 cpassert(csym%kplink(1, i) /= 0)
1157 cpassert(kpop(i) /= 0)
1158 END DO
1159
1160 END SUBROUTINE reduce_kpoint_mesh
1161
1162!> \brief ...
1163!> \param nk ...
1164!> \param xkp ...
1165!> \param wkp ...
1166!> \param shift ...
1167! **************************************************************************************************
1168 SUBROUTINE full_grid_gen(nk, xkp, wkp, shift)
1169 INTEGER, INTENT(IN) :: nk(3)
1170 REAL(kind=dp), DIMENSION(:, :) :: xkp
1171 REAL(kind=dp), DIMENSION(:) :: wkp
1172 REAL(kind=dp), INTENT(IN) :: shift(3)
1173
1174 INTEGER :: i, ix, iy, iz
1175 REAL(kind=dp) :: kpt_latt(3)
1176
1177 wkp = 0.0_dp
1178 i = 0
1179 DO ix = 1, nk(1)
1180 DO iy = 1, nk(2)
1181 DO iz = 1, nk(3)
1182 i = i + 1
1183 kpt_latt(1) = real(2*ix - nk(1) - 1, kind=dp)/(2._dp*real(nk(1), kind=dp))
1184 kpt_latt(2) = real(2*iy - nk(2) - 1, kind=dp)/(2._dp*real(nk(2), kind=dp))
1185 kpt_latt(3) = real(2*iz - nk(3) - 1, kind=dp)/(2._dp*real(nk(3), kind=dp))
1186 xkp(1:3, i) = kpt_latt(1:3)
1187 wkp(i) = 1.0_dp
1188 END DO
1189 END DO
1190 END DO
1191 DO i = 1, nk(1)*nk(2)*nk(3)
1192 xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
1193 END DO
1194
1195 END SUBROUTINE full_grid_gen
1196
1197! **************************************************************************************************
1198!> \brief ...
1199!> \param xkp ...
1200!> \param wkp ...
1201!> \param link ...
1202! **************************************************************************************************
1203 SUBROUTINE inversion_symm(xkp, wkp, link)
1204 REAL(kind=dp), DIMENSION(:, :) :: xkp
1205 REAL(kind=dp), DIMENSION(:) :: wkp
1206 INTEGER, DIMENSION(:) :: link
1207
1208 INTEGER :: i, j, nkpts
1209 REAL(kind=dp), DIMENSION(3) :: diff
1210
1211 nkpts = SIZE(wkp, 1)
1212
1213 link(:) = 0
1214 DO i = 1, nkpts
1215 IF (link(i) == 0) link(i) = i
1216 DO j = i + 1, nkpts
1217 IF (wkp(j) == 0) cycle
1218 diff(1:3) = xkp(1:3, i) + xkp(1:3, j)
1219 diff(1:3) = diff(1:3) - anint(diff(1:3))
1220 IF (all(abs(diff(1:3)) < 1.e-12_dp)) THEN
1221 wkp(i) = wkp(i) + wkp(j)
1222 wkp(j) = 0.0_dp
1223 link(j) = i
1224 EXIT
1225 END IF
1226 END DO
1227 END DO
1228
1229 END SUBROUTINE inversion_symm
1230
1231! **************************************************************************************************
1232!> \brief ...
1233!> \param x ...
1234!> \param r ...
1235!> \return ...
1236! **************************************************************************************************
1237 FUNCTION kp_apply_operation(x, r) RESULT(y)
1238 REAL(kind=dp), INTENT(IN) :: x(3), r(3, 3)
1239 REAL(kind=dp) :: y(3)
1240
1241 y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
1242 y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
1243 y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
1244
1245 END FUNCTION kp_apply_operation
1246
1247! **************************************************************************************************
1248!> \brief ...
1249!> \param csym ...
1250! **************************************************************************************************
1251 SUBROUTINE print_crys_symmetry(csym)
1252 TYPE(csym_type) :: csym
1253
1254 INTEGER :: i, iunit, j, plevel
1255
1256 iunit = csym%punit
1257 IF (iunit >= 0) THEN
1258 plevel = csym%plevel
1259 WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
1260 IF (csym%symlib) THEN
1261 WRITE (iunit, '(A,T71,A10)') " International Symbol: ", adjustr(trim(csym%international_symbol))
1262 WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", adjustr(trim(csym%pointgroup_symbol))
1263 WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", adjustr(trim(csym%schoenflies))
1264 !
1265 WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
1266 IF (plevel > 0) THEN
1267 DO i = 1, csym%n_operations
1268 WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
1269 " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
1270 WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
1271 END DO
1272 END IF
1273 ELSE
1274 WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not availale"
1275 END IF
1276 END IF
1277
1278 END SUBROUTINE print_crys_symmetry
1279
1280! **************************************************************************************************
1281!> \brief ...
1282!> \param csym ...
1283! **************************************************************************************************
1284 SUBROUTINE print_kp_symmetry(csym)
1285 TYPE(csym_type), INTENT(IN) :: csym
1286
1287 INTEGER :: i, iunit, nat, plevel
1288
1289 iunit = csym%punit
1290 IF (iunit >= 0) THEN
1291 plevel = csym%plevel
1292 WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
1293 WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
1294 WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
1295 DO i = 1, csym%nkpoint
1296 WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), nint(csym%wkpoint(i))
1297 END DO
1298 WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
1299 WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
1300 DO i = 1, csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
1301 WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
1302 csym%kplink(1:2, i), csym%kpop(i)
1303 END DO
1304 IF (csym%nrtot > 0) THEN
1305 WRITE (iunit, '(/,A)') " Atom Transformation Table"
1306 nat = SIZE(csym%f0, 1)
1307 DO i = 1, csym%nrtot
1308 WRITE (iunit, '(T10,A,I5,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
1309 END DO
1310 END IF
1311 END IF
1312
1313 END SUBROUTINE print_kp_symmetry
1314
1315END 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::kpoint_gen
subroutine, public kpoint_gen(csym, nk, symm, shift, full_grid, inversion_symmetry_only, use_spglib_reduction, use_spglib_backend)
...
Definition cryssym.F:242
cryssym::print_crys_symmetry
subroutine, public print_crys_symmetry(csym)
...
Definition cryssym.F:1252
cryssym::crys_sym_gen
subroutine, public crys_sym_gen(csym, scoor, types, hmat, delta, iounit)
...
Definition cryssym.F:143
cryssym::release_csym_type
subroutine, public release_csym_type(csym)
Release the CSYM type.
Definition cryssym.F:89
cryssym::print_kp_symmetry
subroutine, public print_kp_symmetry(csym)
...
Definition cryssym.F:1285
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