Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : cite_reference
16 : USE kinds, ONLY: dp
17 : USE kpsym, ONLY: group1s,&
18 : k290s
19 : USE spglib_f08, ONLY: spg_get_international,&
20 : spg_get_major_version,&
21 : spg_get_micro_version,&
22 : spg_get_minor_version,&
23 : spg_get_multiplicity,&
24 : spg_get_pointgroup,&
25 : spg_get_schoenflies,&
26 : spg_get_symmetry
27 : USE string_utilities, ONLY: strip_control_codes
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 : PUBLIC :: csym_type, release_csym_type, print_crys_symmetry, print_kp_symmetry
33 : PUBLIC :: crys_sym_gen, kpoint_gen
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cryssym'
36 :
37 : ! **************************************************************************************************
38 : !> \brief CSM type
39 : !> \par Content:
40 : !>
41 : ! **************************************************************************************************
42 : TYPE csym_type
43 : LOGICAL :: symlib = .FALSE.
44 : LOGICAL :: fullgrid = .FALSE.
45 : INTEGER :: plevel = 0
46 : INTEGER :: punit = -1
47 : INTEGER :: istriz = -1
48 : REAL(KIND=dp) :: delta = 1.0e-8_dp
49 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp
50 : ! KPOINTS
51 : REAL(KIND=dp), DIMENSION(3) :: wvk0 = 0.0_dp
52 : INTEGER, DIMENSION(3) :: mesh = 0
53 : INTEGER :: nkpoint = 0
54 : INTEGER :: nat = 0
55 : INTEGER, DIMENSION(:), ALLOCATABLE :: atype
56 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: scoord
57 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: xkpoint
58 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: wkpoint
59 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: kpmesh
60 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: kplink
61 : INTEGER, DIMENSION(:), ALLOCATABLE :: kpop
62 : !SPGLIB
63 : CHARACTER(len=11) :: international_symbol = ""
64 : CHARACTER(len=6) :: pointgroup_symbol = ""
65 : CHARACTER(len=10) :: schoenflies = ""
66 : INTEGER :: n_operations = 0
67 : INTEGER, DIMENSION(:, :, :), ALLOCATABLE :: rotations
68 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: translations
69 : !K290
70 : REAL(KIND=dp), DIMENSION(3, 3, 48) :: rt = 0.0_dp
71 : REAL(KIND=dp), DIMENSION(3, 48) :: vt = 0.0_dp
72 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0
73 : INTEGER :: nrtot = 0
74 : INTEGER, DIMENSION(48) :: ibrot = 1
75 : END TYPE csym_type
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief Release the CSYM type
81 : !> \param csym The CSYM type
82 : ! **************************************************************************************************
83 445 : SUBROUTINE release_csym_type(csym)
84 : TYPE(csym_type) :: csym
85 :
86 445 : IF (ALLOCATED(csym%rotations)) THEN
87 444 : DEALLOCATE (csym%rotations)
88 : END IF
89 445 : IF (ALLOCATED(csym%translations)) THEN
90 444 : DEALLOCATE (csym%translations)
91 : END IF
92 445 : IF (ALLOCATED(csym%atype)) THEN
93 445 : DEALLOCATE (csym%atype)
94 : END IF
95 445 : IF (ALLOCATED(csym%scoord)) THEN
96 445 : DEALLOCATE (csym%scoord)
97 : END IF
98 445 : IF (ALLOCATED(csym%xkpoint)) THEN
99 152 : DEALLOCATE (csym%xkpoint)
100 : END IF
101 445 : IF (ALLOCATED(csym%wkpoint)) THEN
102 152 : DEALLOCATE (csym%wkpoint)
103 : END IF
104 445 : IF (ALLOCATED(csym%kpmesh)) THEN
105 152 : DEALLOCATE (csym%kpmesh)
106 : END IF
107 445 : IF (ALLOCATED(csym%kplink)) THEN
108 152 : DEALLOCATE (csym%kplink)
109 : END IF
110 445 : IF (ALLOCATED(csym%kpop)) THEN
111 152 : DEALLOCATE (csym%kpop)
112 : END IF
113 445 : IF (ALLOCATED(csym%f0)) THEN
114 0 : DEALLOCATE (csym%f0)
115 : END IF
116 :
117 445 : END SUBROUTINE release_csym_type
118 :
119 : ! **************************************************************************************************
120 : !> \brief ...
121 : !> \param csym ...
122 : !> \param scoor ...
123 : !> \param types ...
124 : !> \param hmat ...
125 : !> \param delta ...
126 : !> \param iounit ...
127 : ! **************************************************************************************************
128 445 : SUBROUTINE crys_sym_gen(csym, scoor, types, hmat, delta, iounit)
129 : TYPE(csym_type) :: csym
130 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scoor
131 : INTEGER, DIMENSION(:), INTENT(IN) :: types
132 : REAL(KIND=dp), INTENT(IN) :: hmat(3, 3)
133 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: delta
134 : INTEGER, INTENT(IN), OPTIONAL :: iounit
135 :
136 : CHARACTER(LEN=*), PARAMETER :: routineN = 'crys_sym_gen'
137 :
138 : INTEGER :: handle, ierr, major, micro, minor, nat, &
139 : nop, tra_mat(3, 3)
140 : LOGICAL :: spglib
141 :
142 445 : CALL timeset(routineN, handle)
143 :
144 : !..total number of atoms
145 445 : nat = SIZE(scoor, 2)
146 445 : csym%nat = nat
147 :
148 : ! output unit
149 445 : IF (PRESENT(iounit)) THEN
150 445 : csym%punit = iounit
151 : ELSE
152 0 : csym%punit = -1
153 : END IF
154 :
155 : ! accuracy for symmetry
156 445 : IF (PRESENT(delta)) THEN
157 445 : csym%delta = delta
158 : ELSE
159 0 : csym%delta = 1.e-6_dp
160 : END IF
161 :
162 : !..set cell values
163 5785 : csym%hmat = hmat
164 :
165 : ! atom types
166 1335 : ALLOCATE (csym%atype(nat))
167 5870 : csym%atype(1:nat) = types(1:nat)
168 :
169 : ! scaled coordinates
170 1335 : ALLOCATE (csym%scoord(3, nat))
171 22145 : csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
172 :
173 445 : csym%n_operations = 0
174 :
175 : !..try spglib
176 445 : major = spg_get_major_version()
177 445 : minor = spg_get_minor_version()
178 445 : micro = spg_get_micro_version()
179 445 : IF (major == 0) THEN
180 0 : CALL cp_warn(__LOCATION__, "Symmetry library SPGLIB not available")
181 0 : spglib = .FALSE.
182 : ELSE
183 445 : spglib = .TRUE.
184 445 : CALL cite_reference(Togo2018)
185 445 : ierr = spg_get_international(csym%international_symbol, TRANSPOSE(hmat), scoor, types, nat, delta)
186 445 : IF (ierr == 0) THEN
187 1 : CALL cp_warn(__LOCATION__, "Symmetry Library SPGLIB failed")
188 1 : spglib = .FALSE.
189 : ELSE
190 444 : nop = spg_get_multiplicity(TRANSPOSE(hmat), scoor, types, nat, delta)
191 2220 : ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
192 444 : csym%n_operations = nop
193 : ierr = spg_get_symmetry(csym%rotations, csym%translations, nop, &
194 444 : TRANSPOSE(hmat), scoor, types, nat, delta)
195 : ! Schoenflies Symbol
196 444 : csym%schoenflies = ' '
197 444 : ierr = spg_get_schoenflies(csym%schoenflies, TRANSPOSE(hmat), scoor, types, nat, delta)
198 : ! Point Group
199 444 : csym%pointgroup_symbol = ' '
200 444 : tra_mat = 0
201 : ierr = spg_get_pointgroup(csym%pointgroup_symbol, tra_mat, &
202 444 : csym%rotations, csym%n_operations)
203 :
204 444 : CALL strip_control_codes(csym%international_symbol)
205 444 : CALL strip_control_codes(csym%schoenflies)
206 444 : CALL strip_control_codes(csym%pointgroup_symbol)
207 : END IF
208 : END IF
209 445 : csym%symlib = spglib
210 :
211 445 : CALL timestop(handle)
212 :
213 445 : END SUBROUTINE crys_sym_gen
214 :
215 : ! **************************************************************************************************
216 : !> \brief ...
217 : !> \param csym ...
218 : !> \param nk ...
219 : !> \param symm ...
220 : !> \param shift ...
221 : !> \param full_grid ...
222 : ! **************************************************************************************************
223 152 : SUBROUTINE kpoint_gen(csym, nk, symm, shift, full_grid)
224 : TYPE(csym_type) :: csym
225 : INTEGER, INTENT(IN) :: nk(3)
226 : LOGICAL, INTENT(IN), OPTIONAL :: symm
227 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: shift(3)
228 : LOGICAL, INTENT(IN), OPTIONAL :: full_grid
229 :
230 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_gen'
231 :
232 : INTEGER :: handle, i, ik, j, nkp, nkpts
233 152 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kpop, xptr
234 : LOGICAL :: fullmesh
235 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkp
236 152 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp
237 :
238 152 : CALL timeset(routineN, handle)
239 :
240 152 : IF (PRESENT(shift)) THEN
241 608 : csym%wvk0 = shift
242 : ELSE
243 0 : csym%wvk0 = 0.0_dp
244 : END IF
245 :
246 152 : csym%istriz = -1
247 152 : IF (PRESENT(symm)) THEN
248 152 : IF (symm) csym%istriz = 1
249 : END IF
250 :
251 152 : IF (PRESENT(full_grid)) THEN
252 152 : fullmesh = full_grid
253 : ELSE
254 : fullmesh = .FALSE.
255 : END IF
256 152 : csym%fullgrid = fullmesh
257 :
258 152 : csym%nkpoint = 0
259 608 : csym%mesh(1:3) = nk(1:3)
260 :
261 152 : nkpts = nk(1)*nk(2)*nk(3)
262 1064 : ALLOCATE (xkp(3, nkpts), wkp(nkpts), kpop(nkpts))
263 : ! kp: link
264 456 : ALLOCATE (csym%kplink(2, nkpts))
265 9944 : csym%kplink = 0
266 :
267 : ! go through all the options
268 152 : IF (csym%symlib) THEN
269 : ! symmetry library is available
270 152 : IF (fullmesh) THEN
271 : ! full mesh requested
272 64 : CALL full_grid_gen(nk, xkp, wkp, shift)
273 64 : IF (csym%istriz == 1) THEN
274 : ! use inversion symmetry
275 48 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
276 : ELSE
277 : ! full kpoint mesh is used
278 : END IF
279 88 : ELSE IF (csym%istriz /= 1) THEN
280 : ! use inversion symmetry
281 88 : CALL full_grid_gen(nk, xkp, wkp, shift)
282 88 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
283 : ELSE
284 : ! use symmetry library to reduce k-points
285 0 : IF (SUM(ABS(csym%wvk0)) /= 0.0_dp) THEN
286 : CALL cp_abort(__LOCATION__, "MacDonald shifted k-point meshes are only "// &
287 0 : "possible without symmetrization.")
288 : END IF
289 :
290 0 : CALL full_grid_gen(nk, xkp, wkp, shift)
291 0 : CALL kp_symmetry(csym, xkp, wkp, kpop)
292 :
293 : END IF
294 : ELSE
295 : ! no symmetry library is available
296 0 : CALL full_grid_gen(nk, xkp, wkp, shift)
297 0 : IF (csym%istriz /= 1 .AND. fullmesh) THEN
298 : ! full kpoint mesh is used
299 0 : DO i = 1, nkpts
300 0 : csym%kplink(1, i) = i
301 : END DO
302 : ELSE
303 : ! use inversion symmetry
304 0 : CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
305 : END IF
306 : END IF
307 : ! count kpoints
308 152 : nkp = 0
309 3416 : DO i = 1, nkpts
310 3416 : IF (wkp(i) > 0.0_dp) nkp = nkp + 1
311 : END DO
312 :
313 : ! store reduced kpoint set
314 152 : csym%nkpoint = nkp
315 760 : ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
316 456 : ALLOCATE (xptr(nkp))
317 3416 : j = 0
318 3416 : DO ik = 1, nkpts
319 3416 : IF (wkp(ik) > 0.0_dp) THEN
320 1700 : j = j + 1
321 1700 : csym%wkpoint(j) = wkp(ik)
322 6800 : csym%xkpoint(1:3, j) = xkp(1:3, ik)
323 1700 : xptr(j) = ik
324 : END IF
325 : END DO
326 152 : CPASSERT(j == nkp)
327 :
328 : ! kp: mesh
329 304 : ALLOCATE (csym%kpmesh(3, nkpts))
330 13208 : csym%kpmesh(1:3, 1:nkpts) = xkp(1:3, 1:nkpts)
331 :
332 : ! kp: link
333 3416 : DO ik = 1, nkpts
334 3264 : i = csym%kplink(1, ik)
335 75342 : DO j = 1, nkp
336 75190 : IF (i == xptr(j)) THEN
337 3164 : csym%kplink(2, ik) = j
338 3164 : EXIT
339 : END IF
340 : END DO
341 : END DO
342 152 : DEALLOCATE (xptr)
343 :
344 : ! kp: operations
345 304 : ALLOCATE (csym%kpop(nkpts))
346 152 : IF (csym%symlib .AND. csym%istriz == 1 .AND. .NOT. fullmesh) THEN
347 : ! atomic symmetry operations possible
348 0 : csym%kpop(1:nkpts) = kpop(1:nkpts)
349 0 : DO ik = 1, nkpts
350 0 : CPASSERT(csym%kpop(ik) /= 0)
351 : END DO
352 : ELSE
353 : ! only time reversal symmetry
354 3416 : DO ik = 1, nkpts
355 3416 : IF (wkp(ik) > 0.0_dp) THEN
356 1700 : csym%kpop(ik) = 1
357 : ELSE
358 1564 : csym%kpop(ik) = 2
359 : END IF
360 : END DO
361 : END IF
362 :
363 152 : DEALLOCATE (xkp, wkp, kpop)
364 :
365 152 : CALL timestop(handle)
366 :
367 152 : END SUBROUTINE kpoint_gen
368 :
369 : ! **************************************************************************************************
370 : !> \brief ...
371 : !> \param csym ...
372 : !> \param xkp ...
373 : !> \param wkp ...
374 : !> \param kpop ...
375 : ! **************************************************************************************************
376 0 : SUBROUTINE kp_symmetry(csym, xkp, wkp, kpop)
377 : TYPE(csym_type) :: csym
378 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
379 : REAL(KIND=dp), DIMENSION(:) :: wkp
380 : INTEGER, DIMENSION(:) :: kpop
381 :
382 : INTEGER :: i, ihc, ihg, ik, indpg, iou, iq1, iq2, &
383 : iq3, istriz, isy, itoj, j, kr, li, lr, &
384 : nat, nc, nhash, nkm, nkp, nkpoint, &
385 : nsp, ntvec
386 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: includ, isc, list, lwght, ty
387 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: f0, lrot
388 : INTEGER, DIMENSION(48) :: ib
389 : REAL(KIND=dp) :: alat
390 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rlist, rx, tvec, wvkl, xkapa
391 : REAL(KIND=dp), DIMENSION(3) :: a1, a2, a3, b1, b2, b3, origin, rr, wvk0
392 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat, strain
393 : REAL(KIND=dp), DIMENSION(3, 3, 48) :: r
394 : REAL(KIND=dp), DIMENSION(3, 48) :: vt
395 :
396 0 : iou = csym%punit
397 0 : hmat = csym%hmat
398 0 : nat = csym%nat
399 0 : iq1 = csym%mesh(1)
400 0 : iq2 = csym%mesh(2)
401 0 : iq3 = csym%mesh(3)
402 : nkpoint = 10*iq1*iq2*iq3
403 0 : nkpoint = 2*MAX(iq1, iq2, iq3)**3
404 0 : wvk0 = csym%wvk0
405 0 : istriz = csym%istriz
406 0 : a1(1:3) = hmat(1:3, 1)
407 0 : a2(1:3) = hmat(1:3, 2)
408 0 : a3(1:3) = hmat(1:3, 3)
409 0 : alat = hmat(1, 1)
410 0 : strain = 0.0_dp
411 0 : ALLOCATE (xkapa(3, nat), rx(3, nat), tvec(3, 200), ty(nat), isc(nat), f0(49, nat))
412 0 : ty(1:nat) = csym%atype(1:nat)
413 0 : nsp = MAXVAL(ty)
414 0 : DO i = 1, nat
415 0 : xkapa(1:3, i) = MATMUL(hmat, csym%scoord(1:3, i))
416 : END DO
417 0 : nhash = 1000
418 0 : ALLOCATE (wvkl(3, nkpoint), rlist(3, nkpoint), includ(nkpoint), list(nhash + nkpoint))
419 0 : ALLOCATE (lrot(48, nkpoint), lwght(nkpoint))
420 :
421 0 : IF (iou > 0) THEN
422 : WRITE (iou, '(/,(T2,A79))') &
423 0 : "*******************************************************************************", &
424 0 : "** Special K-Point Generation by K290 **", &
425 0 : "*******************************************************************************"
426 : END IF
427 :
428 : CALL K290s(iou, nat, nkpoint, nsp, iq1, iq2, iq3, istriz, &
429 : a1, a2, a3, alat, strain, xkapa, rx, tvec, &
430 : ty, isc, f0, ntvec, wvk0, wvkl, lwght, lrot, &
431 0 : nhash, includ, list, rlist, csym%delta)
432 :
433 : CALL GROUP1s(0, a1, a2, a3, nat, ty, xkapa, b1, b2, b3, &
434 : ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
435 0 : vt, f0, r, tvec, origin, rx, isc, csym%delta)
436 :
437 0 : IF (iou > 0) THEN
438 : WRITE (iou, '((T2,A79))') &
439 0 : "*******************************************************************************", &
440 0 : "** Finished K290 **", &
441 0 : "*******************************************************************************"
442 : END IF
443 :
444 0 : csym%rt = r
445 0 : csym%vt = vt
446 0 : csym%nrtot = nc
447 0 : ALLOCATE (csym%f0(nat, nc))
448 0 : DO i = 1, nc
449 0 : csym%f0(1:nat, i) = f0(i, 1:nat)
450 : END DO
451 0 : csym%ibrot = 0
452 0 : csym%ibrot(1:nc) = ib(1:nc)
453 :
454 0 : kpop = 0
455 0 : nkm = iq1*iq2*iq3
456 0 : nkp = 0
457 0 : DO i = 1, nkm
458 0 : IF (lwght(i) == 0) EXIT
459 0 : nkp = nkp + 1
460 : END DO
461 0 : wkp = 0
462 : ik = 0
463 0 : DO i = 1, nkp
464 0 : DO j = 1, nkm
465 0 : wvk0(1:3) = xkp(1:3, j) - wvkl(1:3, i)
466 0 : IF (ALL(ABS(wvk0(1:3)) < 1.e-12_dp)) THEN
467 0 : wkp(j) = lwght(i)
468 0 : itoj = j
469 0 : EXIT
470 : END IF
471 : END DO
472 0 : DO lr = 1, lwght(i)
473 0 : kr = lrot(lr, i)
474 0 : rr(1:3) = kp_apply_operation(wvkl(1:3, i), r(1:3, 1:3, ABS(kr)))
475 0 : IF (kr < 0) rr(1:3) = -rr(1:3)
476 0 : DO j = 1, nkm
477 0 : wvk0(1:3) = xkp(1:3, j) - rr(1:3)
478 0 : IF (ALL(ABS(wvk0(1:3)) < 1.e-12_dp)) THEN
479 0 : csym%kplink(1, j) = itoj
480 0 : kpop(j) = kr
481 0 : EXIT
482 : END IF
483 : END DO
484 : END DO
485 : END DO
486 0 : DEALLOCATE (xkapa, rx, tvec, ty, isc, f0)
487 0 : DEALLOCATE (wvkl, rlist, includ, list)
488 0 : DEALLOCATE (lrot, lwght)
489 :
490 0 : END SUBROUTINE kp_symmetry
491 : ! **************************************************************************************************
492 : !> \brief ...
493 : !> \param nk ...
494 : !> \param xkp ...
495 : !> \param wkp ...
496 : !> \param shift ...
497 : ! **************************************************************************************************
498 152 : SUBROUTINE full_grid_gen(nk, xkp, wkp, shift)
499 : INTEGER, INTENT(IN) :: nk(3)
500 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
501 : REAL(KIND=dp), DIMENSION(:) :: wkp
502 : REAL(KIND=dp), INTENT(IN) :: shift(3)
503 :
504 : INTEGER :: i, ix, iy, iz
505 : REAL(KIND=dp) :: kpt_latt(3)
506 :
507 3416 : wkp = 0.0_dp
508 152 : i = 0
509 526 : DO ix = 1, nk(1)
510 1568 : DO iy = 1, nk(2)
511 4680 : DO iz = 1, nk(3)
512 3264 : i = i + 1
513 3264 : kpt_latt(1) = REAL(2*ix - nk(1) - 1, KIND=dp)/(2._dp*REAL(nk(1), KIND=dp))
514 3264 : kpt_latt(2) = REAL(2*iy - nk(2) - 1, KIND=dp)/(2._dp*REAL(nk(2), KIND=dp))
515 3264 : kpt_latt(3) = REAL(2*iz - nk(3) - 1, KIND=dp)/(2._dp*REAL(nk(3), KIND=dp))
516 13056 : xkp(1:3, i) = kpt_latt(1:3)
517 4306 : wkp(i) = 1.0_dp
518 : END DO
519 : END DO
520 : END DO
521 3416 : DO i = 1, nk(1)*nk(2)*nk(3)
522 13208 : xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
523 : END DO
524 :
525 152 : END SUBROUTINE full_grid_gen
526 :
527 : ! **************************************************************************************************
528 : !> \brief ...
529 : !> \param xkp ...
530 : !> \param wkp ...
531 : !> \param link ...
532 : ! **************************************************************************************************
533 136 : SUBROUTINE inversion_symm(xkp, wkp, link)
534 : REAL(KIND=dp), DIMENSION(:, :) :: xkp
535 : REAL(KIND=dp), DIMENSION(:) :: wkp
536 : INTEGER, DIMENSION(:) :: link
537 :
538 : INTEGER :: i, j, nkpts
539 :
540 136 : nkpts = SIZE(wkp, 1)
541 :
542 3300 : link(:) = 0
543 3300 : DO i = 1, nkpts
544 3164 : IF (link(i) == 0) link(i) = i
545 110072 : DO j = i + 1, nkpts
546 108336 : IF (wkp(j) == 0) CYCLE
547 94432 : IF (ALL(xkp(:, i) == -xkp(:, j))) THEN
548 1564 : wkp(i) = wkp(i) + wkp(j)
549 1564 : wkp(j) = 0.0_dp
550 1564 : link(j) = i
551 1564 : EXIT
552 : END IF
553 : END DO
554 : END DO
555 :
556 136 : END SUBROUTINE inversion_symm
557 :
558 : ! **************************************************************************************************
559 : !> \brief ...
560 : !> \param x ...
561 : !> \param r ...
562 : !> \return ...
563 : ! **************************************************************************************************
564 0 : FUNCTION kp_apply_operation(x, r) RESULT(y)
565 : REAL(KIND=dp), INTENT(IN) :: x(3), r(3, 3)
566 : REAL(KIND=dp) :: y(3)
567 :
568 0 : y(1) = r(1, 1)*x(1) + r(1, 2)*x(2) + r(1, 3)*x(3)
569 0 : y(2) = r(2, 1)*x(1) + r(2, 2)*x(2) + r(2, 3)*x(3)
570 0 : y(3) = r(3, 1)*x(1) + r(3, 2)*x(2) + r(3, 3)*x(3)
571 :
572 0 : END FUNCTION kp_apply_operation
573 :
574 : ! **************************************************************************************************
575 : !> \brief ...
576 : !> \param csym ...
577 : ! **************************************************************************************************
578 341 : SUBROUTINE print_crys_symmetry(csym)
579 : TYPE(csym_type) :: csym
580 :
581 : INTEGER :: i, iunit, j, plevel
582 :
583 341 : iunit = csym%punit
584 341 : IF (iunit >= 0) THEN
585 296 : plevel = csym%plevel
586 296 : WRITE (iunit, "(/,T2,A)") "Crystal Symmetry Information"
587 296 : IF (csym%symlib) THEN
588 295 : WRITE (iunit, '(A,T71,A10)') " International Symbol: ", ADJUSTR(TRIM(csym%international_symbol))
589 295 : WRITE (iunit, '(A,T71,A10)') " Point Group Symbol: ", ADJUSTR(TRIM(csym%pointgroup_symbol))
590 295 : WRITE (iunit, '(A,T71,A10)') " Schoenflies Symbol: ", ADJUSTR(TRIM(csym%schoenflies))
591 : !
592 295 : WRITE (iunit, '(A,T71,I10)') " Number of Symmetry Operations: ", csym%n_operations
593 295 : IF (plevel > 0) THEN
594 0 : DO i = 1, csym%n_operations
595 : WRITE (iunit, '(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
596 0 : " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
597 0 : WRITE (iunit, '(T36,3F15.7)') csym%translations(:, i)
598 : END DO
599 : END IF
600 : ELSE
601 1 : WRITE (iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not availale"
602 : END IF
603 : END IF
604 :
605 341 : END SUBROUTINE print_crys_symmetry
606 :
607 : ! **************************************************************************************************
608 : !> \brief ...
609 : !> \param csym ...
610 : ! **************************************************************************************************
611 48 : SUBROUTINE print_kp_symmetry(csym)
612 : TYPE(csym_type), INTENT(IN) :: csym
613 :
614 : INTEGER :: i, iunit, nat, plevel
615 :
616 48 : iunit = csym%punit
617 48 : IF (iunit >= 0) THEN
618 3 : plevel = csym%plevel
619 3 : WRITE (iunit, "(/,T2,A)") "K-point Symmetry Information"
620 3 : WRITE (iunit, '(A,T67,I14)') " Number of Special K-points: ", csym%nkpoint
621 3 : WRITE (iunit, '(T19,A,T74,A)') " Wavevector Basis ", " Weight"
622 175 : DO i = 1, csym%nkpoint
623 175 : WRITE (iunit, '(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), NINT(csym%wkpoint(i))
624 : END DO
625 3 : WRITE (iunit, '(/,A,T63,3I6)') " K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
626 3 : WRITE (iunit, '(T19,A,T54,A)') " Wavevector Basis ", " Special Points Rotation"
627 347 : DO i = 1, csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
628 344 : WRITE (iunit, '(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
629 691 : csym%kplink(1:2, i), csym%kpop(i)
630 : END DO
631 3 : IF (csym%nrtot > 0) THEN
632 0 : WRITE (iunit, '(/,A)') " Atom Transformation Table"
633 0 : nat = SIZE(csym%f0, 1)
634 0 : DO i = 1, csym%nrtot
635 0 : WRITE (iunit, '(T10,A,I3,(T21,12I5))') " Rot=", csym%ibrot(i), csym%f0(1:nat, i)
636 : END DO
637 : END IF
638 : END IF
639 :
640 48 : END SUBROUTINE print_kp_symmetry
641 :
642 0 : END MODULE cryssym
|
