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