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