22 #include "./base/base_uses.f90"
29 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'cryssym'
37 LOGICAL :: symlib = .false.
38 LOGICAL :: fullgrid = .false.
41 INTEGER :: istriz = -1
42 REAL(KIND=
dp) :: delta = 1.0e-8_dp
43 REAL(KIND=
dp),
DIMENSION(3, 3) :: hmat = 0.0_dp
45 REAL(KIND=
dp),
DIMENSION(3) :: wvk0 = 0.0_dp
46 INTEGER,
DIMENSION(3) :: mesh = 0
47 INTEGER :: nkpoint = 0
49 INTEGER,
DIMENSION(:),
ALLOCATABLE :: atype
50 REAL(KIND=
dp),
DIMENSION(:, :),
ALLOCATABLE :: scoord
51 REAL(KIND=
dp),
DIMENSION(:, :),
ALLOCATABLE :: xkpoint
52 REAL(KIND=
dp),
DIMENSION(:),
ALLOCATABLE :: wkpoint
53 REAL(KIND=
dp),
DIMENSION(:, :),
ALLOCATABLE :: kpmesh
54 INTEGER,
DIMENSION(:, :),
ALLOCATABLE :: kplink
55 INTEGER,
DIMENSION(:),
ALLOCATABLE :: kpop
57 CHARACTER(len=11) :: international_symbol =
""
58 CHARACTER(len=6) :: pointgroup_symbol =
""
59 CHARACTER(len=10) :: schoenflies =
""
60 INTEGER :: n_operations = 0
61 INTEGER,
DIMENSION(:, :, :),
ALLOCATABLE :: rotations
62 REAL(KIND=
dp),
DIMENSION(:, :),
ALLOCATABLE :: translations
72 TYPE(csym_type) :: csym
74 IF (
ALLOCATED(csym%rotations))
THEN
75 DEALLOCATE (csym%rotations)
77 IF (
ALLOCATED(csym%translations))
THEN
78 DEALLOCATE (csym%translations)
80 IF (
ALLOCATED(csym%atype))
THEN
81 DEALLOCATE (csym%atype)
83 IF (
ALLOCATED(csym%scoord))
THEN
84 DEALLOCATE (csym%scoord)
86 IF (
ALLOCATED(csym%xkpoint))
THEN
87 DEALLOCATE (csym%xkpoint)
89 IF (
ALLOCATED(csym%wkpoint))
THEN
90 DEALLOCATE (csym%wkpoint)
92 IF (
ALLOCATED(csym%kpmesh))
THEN
93 DEALLOCATE (csym%kpmesh)
95 IF (
ALLOCATED(csym%kplink))
THEN
96 DEALLOCATE (csym%kplink)
98 IF (
ALLOCATED(csym%kpop))
THEN
99 DEALLOCATE (csym%kpop)
114 TYPE(csym_type) :: csym
115 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: scoor
116 INTEGER,
DIMENSION(:),
INTENT(IN) :: types
117 REAL(kind=
dp),
INTENT(IN) :: hmat(3, 3)
118 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: delta
119 INTEGER,
INTENT(IN),
OPTIONAL :: iounit
121 CHARACTER(LEN=*),
PARAMETER :: routinen =
'crys_sym_gen'
123 INTEGER :: handle, ierr, major, micro, minor, nat, &
127 CALL timeset(routinen, handle)
134 IF (
PRESENT(iounit))
THEN
141 IF (
PRESENT(delta))
THEN
144 csym%delta = 1.e-6_dp
151 ALLOCATE (csym%atype(nat))
152 csym%atype(1:nat) = types(1:nat)
155 ALLOCATE (csym%scoord(3, nat))
156 csym%scoord(1:3, 1:nat) = scoor(1:3, 1:nat)
158 csym%n_operations = 0
165 CALL cp_warn(__location__,
"Symmetry library SPGLIB not available")
172 CALL cp_warn(__location__,
"Symmetry Library SPGLIB failed")
176 ALLOCATE (csym%rotations(3, 3, nop), csym%translations(3, nop))
177 csym%n_operations = nop
179 transpose(hmat), scoor, types, nat, delta)
181 csym%schoenflies =
' '
184 csym%pointgroup_symbol =
' '
187 csym%rotations, csym%n_operations)
196 CALL timestop(handle)
209 TYPE(csym_type) :: csym
210 INTEGER,
INTENT(IN) :: nk(3)
211 LOGICAL,
INTENT(IN),
OPTIONAL :: symm
212 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: shift(3)
213 LOGICAL,
INTENT(IN),
OPTIONAL :: full_grid
215 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_gen'
217 INTEGER :: handle, i, ik, is_shift(3), &
218 is_time_reversal, j, nkp, nkpts
219 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: iwkp, xptr
220 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: ixkp
222 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: wkp
223 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: xkp
224 REAL(kind=
dp),
DIMENSION(3) :: rxkp
226 CALL timeset(routinen, handle)
228 IF (
PRESENT(shift))
THEN
235 IF (
PRESENT(symm))
THEN
236 IF (symm) csym%istriz = 1
239 IF (
PRESENT(full_grid))
THEN
244 csym%fullgrid = fullmesh
247 csym%mesh(1:3) = nk(1:3)
249 nkpts = nk(1)*nk(2)*nk(3)
250 ALLOCATE (xkp(3, nkpts), wkp(nkpts))
252 ALLOCATE (csym%kplink(2, nkpts))
256 IF (csym%symlib)
THEN
260 CALL full_grid_gen(nk, xkp, wkp, shift)
261 IF (csym%istriz == 1)
THEN
263 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
267 ELSE IF (csym%istriz /= 1)
THEN
269 CALL full_grid_gen(nk, xkp, wkp, shift)
270 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
273 IF (sum(abs(csym%wvk0)) /= 0.0_dp)
THEN
274 CALL cp_abort(__location__,
"MacDonald shifted k-point meshes are only "// &
275 "possible without symmetrization.")
277 ALLOCATE (ixkp(3, nkpts), iwkp(nkpts))
278 is_shift(1:3) = mod(nk(1:3) + 1, 2)
281 transpose(csym%hmat), csym%scoord, csym%atype, &
282 csym%nat, csym%delta)
285 xkp(:, i) = real(is_shift + 2*ixkp(:, i), kind=
dp)/real(2*nk(:), kind=
dp)
287 wkp(j) = wkp(j) + 1.0_dp
289 csym%kplink(1, 1:nkpts) = iwkp(1:nkpts) + 1
290 DEALLOCATE (ixkp, iwkp)
294 CALL full_grid_gen(nk, xkp, wkp, shift)
295 IF (csym%istriz /= 1 .AND. fullmesh)
THEN
298 csym%kplink(1, i) = i
302 CALL inversion_symm(xkp, wkp, csym%kplink(1, :))
308 IF (wkp(i) > 0.0_dp) nkp = nkp + 1
313 ALLOCATE (csym%xkpoint(3, nkp), csym%wkpoint(nkp))
317 IF (wkp(ik) > 0.0_dp)
THEN
319 csym%wkpoint(j) = wkp(ik)
320 csym%xkpoint(1:3, j) = xkp(1:3, ik)
327 ALLOCATE (csym%kpmesh(3, nkpts))
328 csym%kpmesh(1:3, 1:nkpts) = xkp(1:3, 1:nkpts)
332 i = csym%kplink(1, ik)
334 IF (i == xptr(j))
THEN
335 csym%kplink(2, ik) = j
343 ALLOCATE (csym%kpop(nkpts))
344 IF (csym%symlib .AND. csym%istriz == 1 .AND. .NOT. fullmesh)
THEN
347 IF (wkp(ik) > 0.0_dp)
THEN
351 j = csym%kplink(2, ik)
352 DO i = 1, csym%n_operations
353 IF (sum(abs(csym%translations(:, i))) > 1.e-10_dp) cycle
354 rxkp(1:3) = kp_apply_operation(csym%xkpoint(1:3, j), csym%rotations(:, :, i))
355 rxkp(1:3) = abs(xkp(1:3, ik) - rxkp(1:3))
356 rxkp(1:3) = rxkp(1:3) - real(nint(rxkp(1:3)), kind=
dp)
357 IF (all((rxkp(1:3)) < 1.e-12_dp))
THEN
362 cpassert(csym%kpop(ik) /= 0)
368 IF (wkp(ik) > 0.0_dp)
THEN
376 DEALLOCATE (xkp, wkp)
378 CALL timestop(handle)
389 SUBROUTINE full_grid_gen(nk, xkp, wkp, shift)
390 INTEGER,
INTENT(IN) :: nk(3)
391 REAL(kind=
dp),
DIMENSION(:, :) :: xkp
392 REAL(kind=
dp),
DIMENSION(:) :: wkp
393 REAL(kind=
dp),
INTENT(IN) :: shift(3)
395 INTEGER :: i, ix, iy, iz
396 REAL(kind=
dp) :: kpt_latt(3)
404 kpt_latt(1) = real(2*ix - nk(1) - 1, kind=
dp)/(2._dp*real(nk(1), kind=
dp))
405 kpt_latt(2) = real(2*iy - nk(2) - 1, kind=
dp)/(2._dp*real(nk(2), kind=
dp))
406 kpt_latt(3) = real(2*iz - nk(3) - 1, kind=
dp)/(2._dp*real(nk(3), kind=
dp))
407 xkp(1:3, i) = kpt_latt(1:3)
412 DO i = 1, nk(1)*nk(2)*nk(3)
413 xkp(1:3, i) = xkp(1:3, i) + shift(1:3)
416 END SUBROUTINE full_grid_gen
424 SUBROUTINE inversion_symm(xkp, wkp, link)
425 REAL(kind=
dp),
DIMENSION(:, :) :: xkp
426 REAL(kind=
dp),
DIMENSION(:) :: wkp
427 INTEGER,
DIMENSION(:) :: link
429 INTEGER :: i, j, nkpts
435 IF (link(i) == 0) link(i) = i
437 IF (wkp(j) == 0) cycle
438 IF (all(xkp(:, i) == -xkp(:, j)))
THEN
439 wkp(i) = wkp(i) + wkp(j)
447 END SUBROUTINE inversion_symm
455 FUNCTION kp_apply_operation(x, r)
RESULT(y)
456 REAL(kind=
dp),
INTENT(IN) :: x(3)
457 INTEGER,
INTENT(IN) :: r(3, 3)
458 REAL(kind=
dp) :: y(3)
460 y(1) = real(r(1, 1),
dp)*x(1) + real(r(1, 2),
dp)*x(2) + real(r(1, 3),
dp)*x(3)
461 y(2) = real(r(2, 1),
dp)*x(1) + real(r(2, 2),
dp)*x(2) + real(r(2, 3),
dp)*x(3)
462 y(3) = real(r(3, 1),
dp)*x(1) + real(r(3, 2),
dp)*x(2) + real(r(3, 3),
dp)*x(3)
464 END FUNCTION kp_apply_operation
473 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: f0
474 TYPE(csym_type) :: csym
475 INTEGER,
INTENT(IN) :: ir
477 INTEGER :: ia, ib, natom
478 REAL(kind=
dp) :: diff
479 REAL(kind=
dp),
DIMENSION(3) :: rb, ri, ro, tr
480 REAL(kind=
dp),
DIMENSION(3, 3) :: rot
483 rot(1:3, 1:3) = csym%rotations(1:3, 1:3, ir)
484 tr(1:3) = csym%translations(1:3, ir)
488 ri(1:3) = csym%scoord(1:3, ia)
489 ro(1) = real(rot(1, 1),
dp)*ri(1) + real(rot(2, 1),
dp)*ri(2) + real(rot(3, 1),
dp)*ri(3) + tr(1)
490 ro(2) = real(rot(1, 2),
dp)*ri(1) + real(rot(2, 2),
dp)*ri(2) + real(rot(3, 2),
dp)*ri(3) + tr(2)
491 ro(3) = real(rot(1, 3),
dp)*ri(1) + real(rot(2, 3),
dp)*ri(2) + real(rot(3, 3),
dp)*ri(3) + tr(3)
493 rb(1:3) = csym%scoord(1:3, ib)
494 diff = sqrt(sum((ri(:) - rb(:))**2))
495 IF (diff < csym%delta)
THEN
500 cpassert(f0(ia) /= 0)
510 TYPE(csym_type) :: csym
512 INTEGER :: i, iunit, j, plevel
517 WRITE (iunit,
"(/,T2,A)")
"Crystal Symmetry Information"
518 IF (csym%symlib)
THEN
519 WRITE (iunit,
'(A,T71,A10)')
" International Symbol: ", adjustr(trim(csym%international_symbol))
520 WRITE (iunit,
'(A,T71,A10)')
" Point Group Symbol: ", adjustr(trim(csym%pointgroup_symbol))
521 WRITE (iunit,
'(A,T71,A10)')
" Schoenflies Symbol: ", adjustr(trim(csym%schoenflies))
523 WRITE (iunit,
'(A,T71,I10)')
" Number of Symmetry Operations: ", csym%n_operations
525 DO i = 1, csym%n_operations
526 WRITE (iunit,
'(A,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
527 " Rotation #: ", i, (csym%rotations(j, :, i), j=1, 3)
528 WRITE (iunit,
'(T36,3F15.7)') csym%translations(:, i)
532 WRITE (iunit,
"(T2,A)")
"SPGLIB for Crystal Symmetry Information determination is not availale"
543 TYPE(csym_type),
INTENT(IN) :: csym
545 INTEGER :: i, iunit, plevel
550 WRITE (iunit,
"(/,T2,A)")
"K-point Symmetry Information"
551 WRITE (iunit,
'(A,T67,I14)')
" Number of Special K-points: ", csym%nkpoint
552 WRITE (iunit,
'(T19,A,T74,A)')
" Wavevector Basis ",
" Weight"
553 DO i = 1, csym%nkpoint
554 WRITE (iunit,
'(T2,i10,3F10.5,T71,I10)') i, csym%xkpoint(1:3, i), nint(csym%wkpoint(i))
556 WRITE (iunit,
'(/,A,T63,3I6)')
" K-point Mesh: ", csym%mesh(1), csym%mesh(2), csym%mesh(3)
557 WRITE (iunit,
'(T19,A,T54,A)')
" Wavevector Basis ",
" Special Points Rotation"
558 DO i = 1, csym%mesh(1)*csym%mesh(2)*csym%mesh(3)
559 WRITE (iunit,
'(T2,i10,3F10.5,T45,3I12)') i, csym%kpmesh(1:3, i), &
560 csym%kplink(1:2, i), csym%kpop(i)
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public togo2018
K-points and crystal symmetry routines.
subroutine, public print_crys_symmetry(csym)
...
subroutine, public crys_sym_gen(csym, scoor, types, hmat, delta, iounit)
...
subroutine, public release_csym_type(csym)
Release the CSYM type.
subroutine, public kpoint_gen(csym, nk, symm, shift, full_grid)
...
subroutine, public apply_rotation_coord(f0, csym, ir)
...
subroutine, public print_kp_symmetry(csym)
...
Defines the basic variable types.
integer, parameter, public dp
Interface for SPGLIB symmetry routines.
integer function, public spg_get_international(symbol, lattice, position, types, num_atom, symprec)
...
integer function, public spg_get_multiplicity(lattice, position, types, num_atom, symprec)
...
integer function, public spg_get_micro_version()
...
integer function, public spg_get_minor_version()
...
integer function, public spg_get_ir_reciprocal_mesh(grid_point, map, mesh, is_shift, is_time_reversal, lattice, position, types, num_atom, symprec)
...
integer function, public spg_get_major_version()
...
integer function, public spg_get_pointgroup(symbol, trans_mat, rotations, num_rotations)
...
integer function, public spg_get_symmetry(rotation, translation, max_size, lattice, position, types, num_atom, symprec)
...
integer function, public spg_get_schoenflies(symbol, lattice, position, types, num_atom, symprec)
...
Utilities for string manipulations.
elemental subroutine, public strip_control_codes(string)
Strip control codes and extended characters from a string, i.e. replace them with blanks.