LCOV - code coverage report
Current view: top level - src - cryssym.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 175 239 73.2 %
Date: 2024-03-28 07:31:50 Functions: 7 11 63.6 %

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

Generated by: LCOV version 1.15