LCOV - code coverage report
Current view: top level - src/motion - space_groups.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 90.2 % 366 330
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 11 11

            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 Space Group Symmetry Module  (version 1.0, January 16, 2020)
      10              : !> \par History
      11              : !>      Pierre-André Cazade [pcazade] 01.2020 - University of Limerick
      12              : !> \author Pierre-André Cazade (first version)
      13              : ! **************************************************************************************************
      14              : MODULE space_groups
      15              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      16              :    USE bibliography,                    ONLY: Togo2018,&
      17              :                                               cite_reference
      18              :    USE cell_methods,                    ONLY: cell_create,&
      19              :                                               init_cell
      20              :    USE cell_types,                      ONLY: cell_copy,&
      21              :                                               cell_type,&
      22              :                                               real_to_scaled,&
      23              :                                               scaled_to_real
      24              :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      25              :                                               cp_subsys_type
      26              :    USE gopt_f_types,                    ONLY: gopt_f_type
      27              :    USE input_constants,                 ONLY: default_cell_direct_id,&
      28              :                                               default_cell_geo_opt_id,&
      29              :                                               default_cell_md_id,&
      30              :                                               default_cell_method_id,&
      31              :                                               default_minimization_method_id,&
      32              :                                               default_ts_method_id
      33              :    USE input_section_types,             ONLY: section_vals_type,&
      34              :                                               section_vals_val_get
      35              :    USE kinds,                           ONLY: dp
      36              :    USE mathlib,                         ONLY: det_3x3,&
      37              :                                               inv_3x3,&
      38              :                                               jacobi
      39              :    USE particle_list_types,             ONLY: particle_list_type
      40              :    USE physcon,                         ONLY: pascal
      41              :    USE space_groups_types,              ONLY: spgr_type
      42              :    USE spglib_f08,                      ONLY: spg_get_international,&
      43              :                                               spg_get_multiplicity,&
      44              :                                               spg_get_pointgroup,&
      45              :                                               spg_get_schoenflies,&
      46              :                                               spg_get_symmetry
      47              :    USE string_utilities,                ONLY: strlcpy_c2f
      48              : #include "../base/base_uses.f90"
      49              : 
      50              :    IMPLICIT NONE
      51              : 
      52              :    PRIVATE
      53              : 
      54              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'space_groups'
      55              : 
      56              :    PUBLIC :: spgr_create, identify_space_group, spgr_find_equivalent_atoms
      57              :    PUBLIC :: spgr_apply_rotations_coord, spgr_apply_rotations_force, print_spgr
      58              :    PUBLIC :: spgr_apply_rotations_stress, spgr_write_stress_tensor
      59              : 
      60              : CONTAINS
      61              : 
      62              : ! **************************************************************************************************
      63              : !> \brief routine creates the space group structure
      64              : !> \param scoor ...
      65              : !> \param types ...
      66              : !> \param cell ...
      67              : !> \param gopt_env ...
      68              : !> \param eps_symmetry ...
      69              : !> \param pol ...
      70              : !> \param ranges ...
      71              : !> \param nparticle ...
      72              : !> \param n_atom ...
      73              : !> \param n_core ...
      74              : !> \param n_shell ...
      75              : !> \param iunit ...
      76              : !> \param print_atoms ...
      77              : !> \par History
      78              : !>      01.2020 created [pcazade]
      79              : !> \author Pierre-André Cazade (first version)
      80              : ! **************************************************************************************************
      81           14 :    SUBROUTINE spgr_create(scoor, types, cell, gopt_env, eps_symmetry, pol, ranges, &
      82              :                           nparticle, n_atom, n_core, n_shell, iunit, print_atoms)
      83              : 
      84              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: scoor
      85              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: types
      86              :       TYPE(cell_type), INTENT(IN), POINTER               :: cell
      87              :       TYPE(gopt_f_type), INTENT(IN), POINTER             :: gopt_env
      88              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: eps_symmetry
      89              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL  :: pol
      90              :       INTEGER, DIMENSION(:, :), INTENT(IN), OPTIONAL     :: ranges
      91              :       INTEGER, INTENT(IN), OPTIONAL                      :: nparticle, n_atom, n_core, n_shell
      92              :       INTEGER, INTENT(IN)                                :: iunit
      93              :       LOGICAL, INTENT(IN)                                :: print_atoms
      94              : 
      95              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_create'
      96              : #ifdef __SPGLIB
      97              :       CHARACTER(LEN=1000)                                :: buffer
      98              :       INTEGER                                            :: ierr, nchars, nop, tra_mat(3, 3)
      99              : #endif
     100              :       INTEGER                                            :: handle, i, j, n_sr_rep
     101           14 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: tmp_types
     102              :       LOGICAL                                            :: spglib
     103           14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: tmp_coor
     104              :       TYPE(spgr_type), POINTER                           :: spgr
     105              : 
     106           14 :       CALL timeset(routineN, handle)
     107              : 
     108           14 :       spgr => gopt_env%spgr
     109           14 :       CPASSERT(ASSOCIATED(spgr))
     110              : 
     111              :       !..total number of particles (atoms plus shells)
     112           14 :       IF (PRESENT(nparticle)) THEN
     113           14 :          CPASSERT(nparticle == SIZE(scoor, 2))
     114           14 :          spgr%nparticle = nparticle
     115              :       ELSE
     116            0 :          spgr%nparticle = SIZE(scoor, 2)
     117              :       END IF
     118              : 
     119           14 :       IF (PRESENT(n_atom)) THEN
     120           14 :          spgr%n_atom = n_atom
     121            0 :       ELSE IF (PRESENT(n_core)) THEN
     122            0 :          spgr%n_atom = spgr%nparticle - n_core
     123            0 :       ELSE IF (PRESENT(n_shell)) THEN
     124            0 :          spgr%n_atom = spgr%nparticle - n_shell
     125              :       ELSE
     126            0 :          spgr%n_atom = spgr%nparticle
     127              :       END IF
     128              : 
     129           14 :       IF (PRESENT(n_core)) THEN
     130           14 :          spgr%n_core = n_core
     131            0 :       ELSE IF (PRESENT(n_shell)) THEN
     132            0 :          spgr%n_core = n_shell
     133              :       END IF
     134              : 
     135           14 :       IF (PRESENT(n_shell)) THEN
     136           14 :          spgr%n_shell = n_shell
     137            0 :       ELSE IF (PRESENT(n_core)) THEN
     138            0 :          spgr%n_shell = n_core
     139              :       END IF
     140              : 
     141           14 :       IF (.NOT. (spgr%nparticle == (spgr%n_atom + spgr%n_shell))) THEN
     142            0 :          CPABORT("spgr_create: nparticle not equal to natom + nshell.")
     143              :       END IF
     144              : 
     145           14 :       spgr%nparticle_sym = spgr%nparticle
     146           14 :       spgr%n_atom_sym = spgr%n_atom
     147           14 :       spgr%n_core_sym = spgr%n_core
     148           14 :       spgr%n_shell_sym = spgr%n_shell
     149              : 
     150           14 :       spgr%iunit = iunit
     151           14 :       spgr%print_atoms = print_atoms
     152              : 
     153              :       ! accuracy for symmetry
     154           14 :       IF (PRESENT(eps_symmetry)) THEN
     155           14 :          spgr%eps_symmetry = eps_symmetry
     156              :       END IF
     157              : 
     158              :       ! vector to test reduced symmetry
     159           14 :       IF (PRESENT(pol)) THEN
     160           14 :          spgr%pol(1) = pol(1)
     161           14 :          spgr%pol(2) = pol(2)
     162           14 :          spgr%pol(3) = pol(3)
     163              :       END IF
     164              : 
     165           42 :       ALLOCATE (spgr%lat(spgr%nparticle))
     166          170 :       spgr%lat = .TRUE.
     167              : 
     168           14 :       IF (PRESENT(ranges)) THEN
     169            0 :          n_sr_rep = SIZE(ranges, 2)
     170            0 :          DO i = 1, n_sr_rep
     171            0 :             DO j = ranges(1, i), ranges(2, i)
     172            0 :                spgr%lat(j) = .FALSE.
     173            0 :                spgr%nparticle_sym = spgr%nparticle_sym - 1
     174            0 :                IF (j <= spgr%n_atom) THEN
     175            0 :                   spgr%n_atom_sym = spgr%n_atom_sym - 1
     176            0 :                ELSE IF (j > spgr%n_atom .AND. j <= spgr%nparticle) THEN
     177            0 :                   spgr%n_core_sym = spgr%n_core_sym - 1
     178            0 :                   spgr%n_shell_sym = spgr%n_shell_sym - 1
     179              :                ELSE
     180            0 :                   CPABORT("Symmetry exclusion range larger than actual number of particles.")
     181              :                END IF
     182              :             END DO
     183              :          END DO
     184              :       END IF
     185              : 
     186           70 :       ALLOCATE (tmp_coor(3, spgr%n_atom_sym), tmp_types(spgr%n_atom_sym))
     187              : 
     188           14 :       j = 0
     189          122 :       DO i = 1, spgr%n_atom
     190          122 :          IF (spgr%lat(i)) THEN
     191          108 :             j = j + 1
     192          432 :             tmp_coor(:, j) = scoor(:, i)
     193          108 :             tmp_types(j) = types(i)
     194              :          END IF
     195              :       END DO
     196              : 
     197              :       !..set cell values
     198           14 :       NULLIFY (spgr%cell_ref)
     199           14 :       CALL cell_create(spgr%cell_ref)
     200           14 :       CALL cell_copy(cell, spgr%cell_ref, tag="CELL_OPT_REF")
     201           14 :       SELECT CASE (gopt_env%type_id)
     202              :       CASE (default_minimization_method_id, default_ts_method_id)
     203            0 :          CALL init_cell(spgr%cell_ref, hmat=cell%hmat)
     204              :       CASE (default_cell_method_id)
     205           28 :          SELECT CASE (gopt_env%cell_method_id)
     206              :          CASE (default_cell_direct_id)
     207           14 :             CALL init_cell(spgr%cell_ref, hmat=gopt_env%h_ref)
     208              :          CASE (default_cell_geo_opt_id)
     209            0 :             CPABORT("SPACE_GROUP_SYMMETRY should not be invoked during the cell step.")
     210              :          CASE (default_cell_md_id)
     211            0 :             CPABORT("SPACE_GROUP_SYMMETRY is not compatible with md.")
     212              :          CASE DEFAULT
     213           14 :             CPABORT("SPACE_GROUP_SYMMETRY invoked with an unknown optimization method.")
     214              :          END SELECT
     215              :       CASE DEFAULT
     216           14 :          CPABORT("SPACE_GROUP_SYMMETRY is not compatible with md.")
     217              :       END SELECT
     218              : 
     219              :       ! atom types
     220           42 :       ALLOCATE (spgr%atype(spgr%nparticle))
     221          170 :       spgr%atype(1:spgr%nparticle) = types(1:spgr%nparticle)
     222              : 
     223           14 :       spgr%n_operations = 0
     224              : 
     225              : #ifdef __SPGLIB
     226           14 :       spglib = .TRUE.
     227           14 :       CALL cite_reference(Togo2018)
     228              :       spgr%space_group_number = spg_get_international(spgr%international_symbol, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
     229          182 :                                                       spgr%n_atom_sym, eps_symmetry)
     230           14 :       buffer = ''
     231           14 :       nchars = strlcpy_c2f(buffer, spgr%international_symbol)
     232           14 :       spgr%international_symbol = buffer(1:nchars)
     233           14 :       IF (spgr%space_group_number == 0) THEN
     234            0 :          CPABORT("Symmetry Library SPGLIB failed, most likely due a problem with the coordinates.")
     235            0 :          spglib = .FALSE.
     236              :       ELSE
     237              :          nop = spg_get_multiplicity(TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
     238           14 :                                     spgr%n_atom_sym, eps_symmetry)
     239           70 :          ALLOCATE (spgr%rotations(3, 3, nop), spgr%translations(3, nop))
     240           56 :          ALLOCATE (spgr%eqatom(nop, spgr%nparticle))
     241           42 :          ALLOCATE (spgr%lop(nop))
     242           14 :          spgr%n_operations = nop
     243          838 :          spgr%lop = .TRUE.
     244              :          ierr = spg_get_symmetry(spgr%rotations, spgr%translations, nop, TRANSPOSE(cell%hmat), &
     245          182 :                                  tmp_coor, tmp_types, spgr%n_atom_sym, eps_symmetry)
     246              :          ! Schoenflies Symbol
     247              :          ierr = spg_get_schoenflies(spgr%schoenflies, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
     248          182 :                                     spgr%n_atom_sym, eps_symmetry)
     249           14 :          buffer = ''
     250           14 :          nchars = strlcpy_c2f(buffer, spgr%schoenflies)
     251           14 :          spgr%schoenflies = buffer(1:nchars)
     252              : 
     253              :          ! Point Group
     254           14 :          tra_mat = 0
     255              :          ierr = spg_get_pointgroup(spgr%pointgroup_symbol, tra_mat, &
     256           14 :                                    spgr%rotations, spgr%n_operations)
     257           14 :          buffer = ''
     258           14 :          nchars = strlcpy_c2f(buffer, spgr%pointgroup_symbol)
     259           14 :          spgr%pointgroup_symbol = buffer(1:nchars)
     260              :       END IF
     261              : #else
     262              :       CPABORT("Symmetry library SPGLIB not available")
     263              :       spglib = .FALSE.
     264              : #endif
     265           14 :       spgr%symlib = spglib
     266              : 
     267           14 :       DEALLOCATE (tmp_coor, tmp_types)
     268              : 
     269           14 :       CALL timestop(handle)
     270              : 
     271           14 :    END SUBROUTINE spgr_create
     272              : 
     273              : ! **************************************************************************************************
     274              : !> \brief routine indentifies the space group and finds rotation matrices.
     275              : !> \param subsys ...
     276              : !> \param geo_section ...
     277              : !> \param gopt_env ...
     278              : !> \param iunit ...
     279              : !> \par History
     280              : !>      01.2020 created [pcazade]
     281              : !> \author Pierre-André Cazade (first version)
     282              : !> \note  rotation matrices innclude translations and translation symmetry:
     283              : !>        it works with supercells as well.
     284              : ! **************************************************************************************************
     285           14 :    SUBROUTINE identify_space_group(subsys, geo_section, gopt_env, iunit)
     286              : 
     287              :       TYPE(cp_subsys_type), INTENT(IN), POINTER          :: subsys
     288              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: geo_section
     289              :       TYPE(gopt_f_type), INTENT(IN), POINTER             :: gopt_env
     290              :       INTEGER, INTENT(IN)                                :: iunit
     291              : 
     292              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'identify_space_group'
     293              : 
     294              :       INTEGER                                            :: handle, i, k, n_atom, n_core, n_shell, &
     295              :                                                             n_sr_rep, nparticle, shell_index
     296           14 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atype
     297           14 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: ranges
     298           14 :       INTEGER, DIMENSION(:), POINTER                     :: tmp
     299              :       LOGICAL                                            :: print_atoms
     300              :       REAL(KIND=dp)                                      :: eps_symmetry
     301           14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: scoord
     302           14 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pol
     303              :       TYPE(cell_type), POINTER                           :: cell
     304              :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     305              :                                                             shell_particles
     306              :       TYPE(spgr_type), POINTER                           :: spgr
     307              : 
     308           14 :       CALL timeset(routineN, handle)
     309              : 
     310           14 :       n_sr_rep = 0
     311           14 :       nparticle = 0
     312              :       n_atom = 0
     313           14 :       n_core = 0
     314           14 :       n_shell = 0
     315              : 
     316           14 :       NULLIFY (particles)
     317           14 :       NULLIFY (core_particles)
     318           14 :       NULLIFY (shell_particles)
     319              : 
     320              :       NULLIFY (cell)
     321           14 :       cell => subsys%cell
     322           14 :       CPASSERT(ASSOCIATED(cell))
     323              : 
     324           14 :       CALL cp_subsys_get(subsys, particles=particles, shell_particles=shell_particles, core_particles=core_particles)
     325              : 
     326           14 :       CPASSERT(ASSOCIATED(particles))
     327           14 :       n_atom = particles%n_els
     328              :       ! Check if we have other kinds of particles in this subsystem
     329           14 :       IF (ASSOCIATED(shell_particles)) THEN
     330            4 :          n_shell = shell_particles%n_els
     331            4 :          CPASSERT(ASSOCIATED(core_particles))
     332            4 :          n_core = subsys%core_particles%n_els
     333              :          ! The same number of shell and core particles is assumed
     334            4 :          CPASSERT(n_core == n_shell)
     335           10 :       ELSE IF (ASSOCIATED(core_particles)) THEN
     336              :          ! This case should not occur at the moment
     337            0 :          CPABORT("Core particles should not be defined without corresponding shell particles.")
     338              :       ELSE
     339              :          n_core = 0
     340              :          n_shell = 0
     341              :       END IF
     342              : 
     343           14 :       nparticle = n_atom + n_shell
     344           70 :       ALLOCATE (scoord(3, nparticle), atype(nparticle))
     345          122 :       DO i = 1, n_atom
     346          108 :          shell_index = particles%els(i)%shell_index
     347          122 :          IF (shell_index == 0) THEN
     348           60 :             CALL real_to_scaled(scoord(1:3, i), particles%els(i)%r(1:3), cell)
     349           60 :             CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, kind_number=atype(i))
     350              :          ELSE
     351           48 :             CALL real_to_scaled(scoord(1:3, i), core_particles%els(shell_index)%r(1:3), cell)
     352           48 :             CALL get_atomic_kind(atomic_kind=core_particles%els(shell_index)%atomic_kind, kind_number=atype(i))
     353           48 :             k = n_atom + shell_index
     354           48 :             CALL real_to_scaled(scoord(1:3, k), shell_particles%els(shell_index)%r(1:3), cell)
     355           48 :             CALL get_atomic_kind(atomic_kind=shell_particles%els(shell_index)%atomic_kind, kind_number=atype(k))
     356              :          END IF
     357              :       END DO
     358              : 
     359           14 :       CALL section_vals_val_get(geo_section, "SPGR_PRINT_ATOMS", l_val=print_atoms)
     360           14 :       CALL section_vals_val_get(geo_section, "EPS_SYMMETRY", r_val=eps_symmetry)
     361           14 :       CALL section_vals_val_get(geo_section, "SYMM_REDUCTION", r_vals=pol)
     362           14 :       CALL section_vals_val_get(geo_section, "SYMM_EXCLUDE_RANGE", n_rep_val=n_sr_rep)
     363           14 :       IF (n_sr_rep > 0) THEN
     364            0 :          ALLOCATE (ranges(2, n_sr_rep))
     365            0 :          DO i = 1, n_sr_rep
     366            0 :             CALL section_vals_val_get(geo_section, "SYMM_EXCLUDE_RANGE", i_rep_val=i, i_vals=tmp)
     367            0 :             ranges(:, i) = tmp(:)
     368              :          END DO
     369              :          CALL spgr_create(scoord, atype, cell, gopt_env, eps_symmetry=eps_symmetry, pol=pol(1:3), &
     370              :                           ranges=ranges, nparticle=nparticle, n_atom=n_atom, &
     371            0 :                           n_core=n_core, n_shell=n_shell, iunit=iunit, print_atoms=print_atoms)
     372            0 :          DEALLOCATE (ranges)
     373              :       ELSE
     374              :          CALL spgr_create(scoord, atype, cell, gopt_env, eps_symmetry=eps_symmetry, pol=pol(1:3), &
     375              :                           nparticle=nparticle, n_atom=n_atom, &
     376           14 :                           n_core=n_core, n_shell=n_shell, iunit=iunit, print_atoms=print_atoms)
     377              :       END IF
     378              : 
     379              :       NULLIFY (spgr)
     380           14 :       spgr => gopt_env%spgr
     381              : 
     382           14 :       CALL spgr_find_equivalent_atoms(spgr, scoord)
     383           14 :       CALL spgr_reduce_symm(spgr)
     384           14 :       CALL spgr_rotations_subset(spgr)
     385              : 
     386           14 :       DEALLOCATE (scoord, atype)
     387              : 
     388           14 :       CALL timestop(handle)
     389              : 
     390           56 :    END SUBROUTINE identify_space_group
     391              : 
     392              : ! **************************************************************************************************
     393              : !> \brief routine indentifies the equivalent atoms for each rotation matrix.
     394              : !> \param spgr ...
     395              : !> \param scoord ...
     396              : !> \par History
     397              : !>      01.2020 created [pcazade]
     398              : !> \author Pierre-André Cazade (first version)
     399              : ! **************************************************************************************************
     400           14 :    SUBROUTINE spgr_find_equivalent_atoms(spgr, scoord)
     401              : 
     402              :       TYPE(spgr_type), INTENT(INOUT), POINTER            :: spgr
     403              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     404              :          INTENT(IN)                                      :: scoord
     405              : 
     406              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_find_equivalent_atoms'
     407              : 
     408              :       INTEGER                                            :: handle, i, ia, ib, ir, j, natom, nop, &
     409              :                                                             nshell
     410              :       REAL(KIND=dp)                                      :: diff
     411              :       REAL(KIND=dp), DIMENSION(3)                        :: rb, ri, ro, tr
     412              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot
     413              : 
     414           14 :       CALL timeset(routineN, handle)
     415              : 
     416           14 :       nop = spgr%n_operations
     417           14 :       natom = spgr%n_atom
     418           14 :       nshell = spgr%n_shell
     419              : 
     420           14 :       IF (.NOT. (spgr%nparticle == (natom + nshell))) THEN
     421            0 :          CPABORT("spgr_find_equivalent_atoms: nparticle not equal to natom + nshell.")
     422              :       END IF
     423              : 
     424          170 :       DO ia = 1, spgr%nparticle
     425         3906 :          spgr%eqatom(:, ia) = ia
     426              :       END DO
     427              : 
     428           14 :       !$OMP PARALLEL DO PRIVATE (ia,ib,ir,ri,rb,ro,rot,tr,diff) SHARED (spgr,scoord,natom,nop) DEFAULT(NONE)
     429              :       DO ia = 1, natom
     430              :          IF (.NOT. spgr%lat(ia)) CYCLE
     431              :          ri(1:3) = scoord(1:3, ia)
     432              :          DO ir = 1, nop
     433              :             rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
     434              :             tr(1:3) = spgr%translations(1:3, ir)
     435              :             DO ib = 1, natom
     436              :                IF (.NOT. spgr%lat(ib)) CYCLE
     437              :                rb(1:3) = scoord(1:3, ib)
     438              :                ro(1) = REAL(rot(1, 1), dp)*rb(1) + REAL(rot(2, 1), dp)*rb(2) + REAL(rot(3, 1), dp)*rb(3) + tr(1)
     439              :                ro(2) = REAL(rot(1, 2), dp)*rb(1) + REAL(rot(2, 2), dp)*rb(2) + REAL(rot(3, 2), dp)*rb(3) + tr(2)
     440              :                ro(3) = REAL(rot(1, 3), dp)*rb(1) + REAL(rot(2, 3), dp)*rb(2) + REAL(rot(3, 3), dp)*rb(3) + tr(3)
     441              :                ro(1) = ro(1) - REAL(NINT(ro(1) - ri(1)), dp)
     442              :                ro(2) = ro(2) - REAL(NINT(ro(2) - ri(2)), dp)
     443              :                ro(3) = ro(3) - REAL(NINT(ro(3) - ri(3)), dp)
     444              :                diff = NORM2(ri(:) - ro(:))
     445              :                IF ((diff < spgr%eps_symmetry) .AND. (spgr%atype(ia) == spgr%atype(ib))) THEN
     446              :                   spgr%eqatom(ir, ia) = ib
     447              :                   EXIT
     448              :                END IF
     449              :             END DO
     450              :          END DO
     451              :       END DO
     452              :       !$OMP END PARALLEL DO
     453              : 
     454           14 :       !$OMP PARALLEL DO PRIVATE (i,j,ia,ib,ir,ri,rb,ro,rot,tr,diff) SHARED (spgr,scoord,natom,nshell,nop) DEFAULT(NONE)
     455              :       DO i = 1, nshell
     456              :          ia = natom + i
     457              :          IF (.NOT. spgr%lat(ia)) CYCLE
     458              :          ri(1:3) = scoord(1:3, ia)
     459              :          DO ir = 1, nop
     460              :             rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
     461              :             tr(1:3) = spgr%translations(1:3, ir)
     462              :             DO j = 1, nshell
     463              :                ib = natom + j
     464              :                IF (.NOT. spgr%lat(ib)) CYCLE
     465              :                rb(1:3) = scoord(1:3, ib)
     466              :                ro(1) = REAL(rot(1, 1), dp)*rb(1) + REAL(rot(2, 1), dp)*rb(2) + REAL(rot(3, 1), dp)*rb(3) + tr(1)
     467              :                ro(2) = REAL(rot(1, 2), dp)*rb(1) + REAL(rot(2, 2), dp)*rb(2) + REAL(rot(3, 2), dp)*rb(3) + tr(2)
     468              :                ro(3) = REAL(rot(1, 3), dp)*rb(1) + REAL(rot(2, 3), dp)*rb(2) + REAL(rot(3, 3), dp)*rb(3) + tr(3)
     469              :                ro(1) = ro(1) - REAL(NINT(ro(1) - ri(1)), dp)
     470              :                ro(2) = ro(2) - REAL(NINT(ro(2) - ri(2)), dp)
     471              :                ro(3) = ro(3) - REAL(NINT(ro(3) - ri(3)), dp)
     472              :                diff = NORM2(ri(:) - ro(:))
     473              :                IF ((diff < spgr%eps_symmetry) .AND. (spgr%atype(ia) == spgr%atype(ib))) THEN
     474              :                   spgr%eqatom(ir, ia) = ib
     475              :                   EXIT
     476              :                END IF
     477              :             END DO
     478              :          END DO
     479              :       END DO
     480              :       !$OMP END PARALLEL DO
     481              : 
     482           14 :       CALL timestop(handle)
     483              : 
     484           14 :    END SUBROUTINE spgr_find_equivalent_atoms
     485              : 
     486              : ! **************************************************************************************************
     487              : !> \brief routine looks for operations compatible with efield
     488              : !> \param spgr ...
     489              : !> \par History
     490              : !>      01.2020 created [pcazade]
     491              : !> \author Pierre-André Cazade (first version)
     492              : ! **************************************************************************************************
     493           14 :    SUBROUTINE spgr_reduce_symm(spgr)
     494              : 
     495              :       TYPE(spgr_type), INTENT(INOUT), POINTER            :: spgr
     496              : 
     497              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'spgr_reduce_symm'
     498              : 
     499              :       INTEGER                                            :: handle, ia, ib, ir, ja, jb, nop, nops, &
     500              :                                                             nparticle
     501           14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: x, xold
     502              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, ro
     503              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot
     504              : 
     505           14 :       CALL timeset(routineN, handle)
     506              : 
     507           14 :       nop = spgr%n_operations
     508           14 :       nparticle = spgr%nparticle
     509           56 :       ALLOCATE (x(3*nparticle), xold(3*nparticle))
     510          482 :       x = 0.0_dp
     511          170 :       DO ia = 1, nparticle
     512          156 :          ja = 3*(ia - 1)
     513          156 :          x(ja + 1) = x(ja + 1) + spgr%pol(1)
     514          156 :          x(ja + 2) = x(ja + 2) + spgr%pol(2)
     515          170 :          x(ja + 3) = x(ja + 3) + spgr%pol(3)
     516              :       END DO
     517          482 :       xold(:) = x(:)
     518              : 
     519              :       nops = 0
     520          838 :       DO ir = 1, nop
     521        12032 :          x = 0.d0
     522          824 :          spgr%lop(ir) = .TRUE.
     523        10712 :          rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
     524         4560 :          DO ia = 1, nparticle
     525         3736 :             IF (.NOT. spgr%lat(ia)) CYCLE
     526         3736 :             ja = 3*(ia - 1)
     527        14944 :             ri(1:3) = xold(ja + 1:ja + 3)
     528         3736 :             ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3)
     529         3736 :             ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3)
     530         3736 :             ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3)
     531        15768 :             x(ja + 1:ja + 3) = ro(1:3)
     532              :          END DO
     533         4560 :          DO ia = 1, nparticle
     534         3736 :             IF (.NOT. spgr%lat(ia)) CYCLE
     535         3736 :             ib = spgr%eqatom(ir, ia)
     536         3736 :             ja = 3*(ia - 1)
     537         3736 :             jb = 3*(ib - 1)
     538        14944 :             ro = x(jb + 1:jb + 3) - xold(ja + 1:ja + 3)
     539              :             spgr%lop(ir) = (spgr%lop(ir) .AND. (ABS(ro(1)) < spgr%eps_symmetry) &
     540              :                             .AND. (ABS(ro(2)) < spgr%eps_symmetry) &
     541         4560 :                             .AND. (ABS(ro(3)) < spgr%eps_symmetry))
     542              :          END DO
     543          838 :          IF (spgr%lop(ir)) nops = nops + 1
     544              :       END DO
     545              : 
     546           14 :       spgr%n_reduced_operations = nops
     547              : 
     548           14 :       DEALLOCATE (x, xold)
     549           14 :       CALL timestop(handle)
     550              : 
     551           14 :    END SUBROUTINE spgr_reduce_symm
     552              : 
     553              : ! **************************************************************************************************
     554              : !> \brief routine looks for unique rotations
     555              : !> \param spgr ...
     556              : !> \par History
     557              : !>      01.2020 created [pcazade]
     558              : !> \author Pierre-André Cazade (first version)
     559              : ! **************************************************************************************************
     560              : 
     561           14 :    SUBROUTINE spgr_rotations_subset(spgr)
     562              : 
     563              :       TYPE(spgr_type), INTENT(INOUT), POINTER            :: spgr
     564              : 
     565              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_rotations_subset'
     566              : 
     567              :       INTEGER                                            :: handle, i, j
     568              :       INTEGER, DIMENSION(3, 3)                           :: d
     569           14 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: mask
     570              : 
     571           14 :       CALL timeset(routineN, handle)
     572              : 
     573           42 :       ALLOCATE (mask(spgr%n_operations))
     574          838 :       mask = .TRUE.
     575              : 
     576          838 :       DO i = 1, spgr%n_operations
     577          838 :          IF (.NOT. spgr%lop(i)) mask(i) = .FALSE.
     578              :       END DO
     579              : 
     580          824 :       DO i = 1, spgr%n_operations - 1
     581          810 :          IF (.NOT. mask(i)) CYCLE
     582        32552 :          DO j = i + 1, spgr%n_operations
     583        32300 :             IF (.NOT. mask(j)) CYCLE
     584       243932 :             d(:, :) = spgr%rotations(:, :, j) - spgr%rotations(:, :, i)
     585       244742 :             IF (SUM(ABS(d)) == 0) mask(j) = .FALSE.
     586              :          END DO
     587              :       END DO
     588              : 
     589           14 :       spgr%n_operations_subset = 0
     590          838 :       DO i = 1, spgr%n_operations
     591          838 :          IF (mask(i)) spgr%n_operations_subset = spgr%n_operations_subset + 1
     592              :       END DO
     593              : 
     594           42 :       ALLOCATE (spgr%rotations_subset(3, 3, spgr%n_operations_subset))
     595              : 
     596           14 :       j = 0
     597          838 :       DO i = 1, spgr%n_operations
     598          838 :          IF (mask(i)) THEN
     599          248 :             j = j + 1
     600         3224 :             spgr%rotations_subset(:, :, j) = spgr%rotations(:, :, i)
     601              :          END IF
     602              :       END DO
     603              : 
     604           14 :       DEALLOCATE (mask)
     605           14 :       CALL timestop(handle)
     606              : 
     607           14 :    END SUBROUTINE spgr_rotations_subset
     608              : 
     609              : ! **************************************************************************************************
     610              : !> \brief routine applies the rotation matrices to the coordinates.
     611              : !> \param spgr ...
     612              : !> \param coord ...
     613              : !> \par History
     614              : !>      01.2020 created [pcazade]
     615              : !> \author Pierre-André Cazade (first version)
     616              : ! **************************************************************************************************
     617           46 :    SUBROUTINE spgr_apply_rotations_coord(spgr, coord)
     618              : 
     619              :       TYPE(spgr_type), INTENT(IN), POINTER               :: spgr
     620              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: coord
     621              : 
     622              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_apply_rotations_coord'
     623              : 
     624              :       INTEGER                                            :: handle, ia, ib, ir, ja, jb, nop, nops, &
     625              :                                                             nparticle
     626           46 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cold
     627              :       REAL(KIND=dp), DIMENSION(3)                        :: rf, ri, rn, ro, tr
     628              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot
     629              : 
     630           46 :       CALL timeset(routineN, handle)
     631              : 
     632          138 :       ALLOCATE (cold(SIZE(coord)))
     633         2074 :       cold(:) = coord(:)
     634              : 
     635           46 :       nop = spgr%n_operations
     636           46 :       nparticle = spgr%nparticle
     637           46 :       nops = spgr%n_reduced_operations
     638              : 
     639           46 :       !$OMP PARALLEL DO PRIVATE (ia,ib,ja,jb,ir,ri,ro,rf,rn,rot,tr) SHARED (spgr,coord,nparticle,nop,nops) DEFAULT(NONE)
     640              :       DO ia = 1, nparticle
     641              :          IF (.NOT. spgr%lat(ia)) CYCLE
     642              :          ja = 3*(ia - 1)
     643              :          CALL real_to_scaled(rf(1:3), coord(ja + 1:ja + 3), spgr%cell_ref)
     644              :          rn(1:3) = 0.d0
     645              :          DO ir = 1, nop
     646              :             IF (.NOT. spgr%lop(ir)) CYCLE
     647              :             ib = spgr%eqatom(ir, ia)
     648              :             rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
     649              :             tr(1:3) = spgr%translations(1:3, ir)
     650              :             jb = 3*(ib - 1)
     651              :             CALL real_to_scaled(ri(1:3), coord(jb + 1:jb + 3), spgr%cell_ref)
     652              :             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)
     653              :             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)
     654              :             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)
     655              :             ro(1) = ro(1) - REAL(NINT(ro(1) - rf(1)), dp)
     656              :             ro(2) = ro(2) - REAL(NINT(ro(2) - rf(2)), dp)
     657              :             ro(3) = ro(3) - REAL(NINT(ro(3) - rf(3)), dp)
     658              :             rn(1:3) = rn(1:3) + ro(1:3)
     659              :          END DO
     660              :          rn = rn/REAL(nops, dp)
     661              :          CALL scaled_to_real(coord(ja + 1:ja + 3), rn(1:3), spgr%cell_ref)
     662              :       END DO
     663              :       !$OMP END PARALLEL DO
     664              : 
     665           46 :       DEALLOCATE (cold)
     666           46 :       CALL timestop(handle)
     667              : 
     668           46 :    END SUBROUTINE spgr_apply_rotations_coord
     669              : 
     670              : ! **************************************************************************************************
     671              : !> \brief routine applies the rotation matrices to the forces.
     672              : !> \param spgr ...
     673              : !> \param force ...
     674              : !> \par History
     675              : !>      01.2020 created [pcazade]
     676              : !> \author Pierre-André Cazade (first version)
     677              : ! **************************************************************************************************
     678           57 :    SUBROUTINE spgr_apply_rotations_force(spgr, force)
     679              : 
     680              :       TYPE(spgr_type), INTENT(IN), POINTER               :: spgr
     681              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: force
     682              : 
     683              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_apply_rotations_force'
     684              : 
     685              :       INTEGER                                            :: handle, ia, ib, ir, ja, jb, nop, nops, &
     686              :                                                             nparticle
     687           57 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: fold
     688              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rn, ro
     689              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: rot
     690              : 
     691           57 :       CALL timeset(routineN, handle)
     692              : 
     693          171 :       ALLOCATE (fold(SIZE(force)))
     694         2943 :       fold(:) = force(:)
     695              : 
     696           57 :       nop = spgr%n_operations
     697           57 :       nparticle = spgr%nparticle
     698           57 :       nops = spgr%n_reduced_operations
     699              : 
     700           57 :       !$OMP PARALLEL DO PRIVATE (ia,ib,ja,jb,ir,ri,ro,rn,rot) SHARED (spgr,force,nparticle,nop,nops) DEFAULT(NONE)
     701              :       DO ia = 1, nparticle
     702              :          IF (.NOT. spgr%lat(ia)) CYCLE
     703              :          ja = 3*(ia - 1)
     704              :          rn(1:3) = 0.d0
     705              :          DO ir = 1, nop
     706              :             IF (.NOT. spgr%lop(ir)) CYCLE
     707              :             ib = spgr%eqatom(ir, ia)
     708              :             rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
     709              :             jb = 3*(ib - 1)
     710              :             CALL real_to_scaled(ri(1:3), force(jb + 1:jb + 3), spgr%cell_ref)
     711              :             ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3)
     712              :             ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3)
     713              :             ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3)
     714              :             rn(1:3) = rn(1:3) + ro(1:3)
     715              :          END DO
     716              :          rn = rn/REAL(nops, dp)
     717              :          CALL scaled_to_real(force(ja + 1:ja + 3), rn(1:3), spgr%cell_ref)
     718              :       END DO
     719              :       !$OMP END PARALLEL DO
     720              : 
     721           57 :       DEALLOCATE (fold)
     722           57 :       CALL timestop(handle)
     723              : 
     724           57 :    END SUBROUTINE spgr_apply_rotations_force
     725              : 
     726              : ! **************************************************************************************************
     727              : !> \brief ...
     728              : !> \param roti ...
     729              : !> \param roto ...
     730              : !> \param nop ...
     731              : !> \param h1 ...
     732              : !> \param h2 ...
     733              : ! **************************************************************************************************
     734           20 :    SUBROUTINE spgr_change_basis(roti, roto, nop, h1, h2)
     735              : 
     736              :       INTEGER, DIMENSION(:, :, :)                        :: roti
     737              :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: roto
     738              :       INTEGER                                            :: nop
     739              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: h1, h2
     740              : 
     741              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'spgr_change_basis'
     742              : 
     743              :       INTEGER                                            :: handle, ir
     744              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: h1ih2, h2ih1, ih1, ih2, r, s
     745              : 
     746           20 :       CALL timeset(routineN, handle)
     747              : 
     748           20 :       ih1 = inv_3x3(h1)
     749           20 :       ih2 = inv_3x3(h2)
     750          800 :       h2ih1 = MATMUL(h2, ih1)
     751          800 :       h1ih2 = MATMUL(h1, ih2)
     752              : 
     753          196 :       DO ir = 1, nop
     754         2288 :          r(:, :) = roti(:, :, ir)
     755         7040 :          s = MATMUL(h2ih1, r)
     756         7040 :          r = MATMUL(s, h1ih2)
     757         2308 :          roto(:, :, ir) = r(:, :)
     758              :       END DO
     759              : 
     760           20 :       CALL timestop(handle)
     761              : 
     762           20 :    END SUBROUTINE spgr_change_basis
     763              : 
     764              : ! **************************************************************************************************
     765              : !> \brief routine applies the rotation matrices to the stress tensor.
     766              : !> \param spgr ...
     767              : !> \param cell ...
     768              : !> \param stress ...
     769              : !> \par History
     770              : !>      01.2020 created [pcazade]
     771              : !> \author Pierre-André Cazade (first version)
     772              : ! **************************************************************************************************
     773           20 :    SUBROUTINE spgr_apply_rotations_stress(spgr, cell, stress)
     774              : 
     775              :       TYPE(spgr_type), INTENT(IN), POINTER               :: spgr
     776              :       TYPE(cell_type), INTENT(IN), POINTER               :: cell
     777              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: stress
     778              : 
     779              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_apply_rotations_stress'
     780              : 
     781              :       INTEGER                                            :: handle, i, ir, j, k, l, nop
     782           20 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: roto
     783              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat1, hmat2, r, stin
     784              : 
     785           20 :       CALL timeset(routineN, handle)
     786              : 
     787          260 :       hmat1 = TRANSPOSE(cell%hmat)
     788              : 
     789           20 :       hmat2 = 0d0
     790           20 :       hmat2(1, 1) = 1.d0
     791           20 :       hmat2(2, 2) = 1.d0
     792           20 :       hmat2(3, 3) = 1.d0
     793              : 
     794           20 :       nop = spgr%n_operations_subset
     795              : 
     796           60 :       ALLOCATE (roto(3, 3, nop))
     797              : 
     798           20 :       CALL spgr_change_basis(spgr%rotations_subset, roto, spgr%n_operations_subset, hmat1, hmat2)
     799              : 
     800           20 :       stin = stress
     801           20 :       stress = 0.d0
     802          196 :       DO ir = 1, nop
     803         2288 :          r(:, :) = roto(:, :, ir)
     804          724 :          DO i = 1, 3
     805         2288 :             DO j = 1, 3
     806         6864 :                DO k = 1, 3
     807        20592 :                   DO l = 1, 3
     808        19008 :                      stress(i, j) = stress(i, j) + (r(k, i)*r(l, j)*stin(k, l))
     809              :                   END DO
     810              :                END DO
     811              :             END DO
     812              :          END DO
     813              :       END DO
     814          260 :       stress = stress/REAL(nop, dp)
     815              : 
     816           20 :       DEALLOCATE (roto)
     817              : 
     818           20 :       CALL timestop(handle)
     819              : 
     820           20 :    END SUBROUTINE spgr_apply_rotations_stress
     821              : 
     822              : ! **************************************************************************************************
     823              : !> \brief routine prints Space Group Information.
     824              : !> \param spgr ...
     825              : !> \par History
     826              : !>      01.2020 created [pcazade]
     827              : !> \author Pierre-André Cazade (first version)
     828              : ! **************************************************************************************************
     829           14 :    SUBROUTINE print_spgr(spgr)
     830              : 
     831              :       TYPE(spgr_type), INTENT(IN), POINTER               :: spgr
     832              : 
     833              :       INTEGER                                            :: i, j
     834              : 
     835           14 :       IF (spgr%iunit > 0) THEN
     836            7 :          WRITE (spgr%iunit, '(/,T2,A,A)') "----------------------------------------", &
     837           14 :             "---------------------------------------"
     838            7 :          WRITE (spgr%iunit, "(T2,A,T25,A,T77,A)") "----", "SPACE GROUP SYMMETRY INFORMATION", "----"
     839            7 :          WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
     840           14 :             "---------------------------------------"
     841            7 :          IF (spgr%symlib) THEN
     842            7 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| SPACE GROUP NUMBER:", &
     843           14 :                spgr%space_group_number
     844            7 :             WRITE (spgr%iunit, '(T2,A,T70,A11)') "SPGR| INTERNATIONAL SYMBOL:", &
     845           14 :                TRIM(ADJUSTR(spgr%international_symbol))
     846            7 :             WRITE (spgr%iunit, '(T2,A,T75,A6)') "SPGR| POINT GROUP SYMBOL:", &
     847           14 :                TRIM(ADJUSTR(spgr%pointgroup_symbol))
     848            7 :             WRITE (spgr%iunit, '(T2,A,T74,A7)') "SPGR| SCHOENFLIES SYMBOL:", &
     849           14 :                TRIM(ADJUSTR(spgr%schoenflies))
     850            7 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF SYMMETRY OPERATIONS:", &
     851           14 :                spgr%n_operations
     852            7 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF UNIQUE ROTATIONS:", &
     853           14 :                spgr%n_operations_subset
     854            7 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF REDUCED SYMMETRY OPERATIONS:", &
     855           14 :                spgr%n_reduced_operations
     856            7 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF PARTICLES AND SYMMETRIZED PARTICLES:", &
     857           14 :                spgr%nparticle, spgr%nparticle_sym
     858            7 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF ATOMS AND SYMMETRIZED ATOMS:", &
     859           14 :                spgr%n_atom, spgr%n_atom_sym
     860            7 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF CORES AND SYMMETRIZED CORES:", &
     861           14 :                spgr%n_core, spgr%n_core_sym
     862            7 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF SHELLS AND SYMMETRIZED SHELLS:", &
     863           14 :                spgr%n_shell, spgr%n_shell_sym
     864            7 :             IF (spgr%print_atoms) THEN
     865            1 :                WRITE (spgr%iunit, *) "SPGR| ACTIVE REDUCED SYMMETRY OPERATIONS:", spgr%lop
     866            1 :                WRITE (spgr%iunit, '(/,T2,A,A)') "----------------------------------------", &
     867            2 :                   "---------------------------------------"
     868            1 :                WRITE (spgr%iunit, '(T2,A,T34,A,T77,A)') "----", "EQUIVALENT ATOMS", "----"
     869            1 :                WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
     870            2 :                   "---------------------------------------"
     871           25 :                DO i = 1, spgr%nparticle
     872          121 :                   DO j = 1, spgr%n_operations
     873           96 :                      WRITE (spgr%iunit, '(T2,A,T52,I8,I8,I8)') "SPGR| ATOM | SYMMETRY OPERATION | EQUIVALENT ATOM", &
     874          216 :                         i, j, spgr%eqatom(j, i)
     875              :                   END DO
     876              :                END DO
     877            1 :                WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
     878            2 :                   "---------------------------------------"
     879            5 :                DO i = 1, spgr%n_operations
     880              :                   WRITE (spgr%iunit, '(T2,A,T46,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
     881           16 :                      "SPGR| SYMMETRY OPERATION #:", i, (spgr%rotations(j, :, i), j=1, 3)
     882            5 :                   WRITE (spgr%iunit, '(T51,3F10.5)') spgr%translations(:, i)
     883              :                END DO
     884              :             END IF
     885              :          ELSE
     886            0 :             WRITE (spgr%iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not availale"
     887              :          END IF
     888              :       END IF
     889              : 
     890           14 :    END SUBROUTINE print_spgr
     891              : 
     892              : ! **************************************************************************************************
     893              : !> \brief Variable precision output of the symmetrized stress tensor
     894              : !>
     895              : !> \param stress tensor ...
     896              : !> \param spgr ...
     897              : !> \par History
     898              : !>      07.2020 adapted to spgr [pcazade]
     899              : !> \author MK (26.08.2010).
     900              : ! **************************************************************************************************
     901           20 :    SUBROUTINE spgr_write_stress_tensor(stress, spgr)
     902              : 
     903              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: stress
     904              :       TYPE(spgr_type), INTENT(IN), POINTER               :: spgr
     905              : 
     906              :       REAL(KIND=dp), DIMENSION(3)                        :: eigval
     907              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: eigvec, stress_tensor
     908              : 
     909          260 :       stress_tensor(:, :) = stress(:, :)*pascal*1.0E-9_dp
     910              : 
     911           20 :       IF (spgr%iunit > 0) THEN
     912              :          WRITE (UNIT=spgr%iunit, FMT='(/,T2,A)') &
     913           11 :             'SPGR STRESS| Symmetrized stress tensor [GPa]'
     914              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T19,3(19X,A1))') &
     915           11 :             'SPGR STRESS|', 'x', 'y', 'z'
     916              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
     917           11 :             'SPGR STRESS|      x', stress_tensor(1, 1:3)
     918              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
     919           11 :             'SPGR STRESS|      y', stress_tensor(2, 1:3)
     920              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
     921           11 :             'SPGR STRESS|      z', stress_tensor(3, 1:3)
     922              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T66,ES20.11)') &
     923           11 :             'SPGR STRESS| 1/3 Trace', (stress_tensor(1, 1) + &
     924              :                                        stress_tensor(2, 2) + &
     925           22 :                                        stress_tensor(3, 3))/3.0_dp
     926              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T66,ES20.11)') &
     927           11 :             'SPGR STRESS| Determinant', det_3x3(stress_tensor(1:3, 1), &
     928              :                                                 stress_tensor(1:3, 2), &
     929           22 :                                                 stress_tensor(1:3, 3))
     930           11 :          eigval(:) = 0.0_dp
     931           11 :          eigvec(:, :) = 0.0_dp
     932           11 :          CALL jacobi(stress_tensor, eigval, eigvec)
     933              :          WRITE (UNIT=spgr%iunit, FMT='(/,T2,A)') &
     934           11 :             'SPGR STRESS| Eigenvectors and eigenvalues of the symmetrized stress tensor [GPa]'
     935              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T19,3(1X,I19))') &
     936           11 :             'SPGR STRESS|', 1, 2, 3
     937              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
     938           11 :             'SPGR STRESS| Eigenvalues', eigval(1:3)
     939              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
     940           11 :             'SPGR STRESS|      x', eigvec(1, 1:3)
     941              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
     942           11 :             'SPGR STRESS|      y', eigvec(2, 1:3)
     943              :          WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
     944           11 :             'SPGR STRESS|      z', eigvec(3, 1:3)
     945              :       END IF
     946              : 
     947           20 :    END SUBROUTINE spgr_write_stress_tensor
     948              : 
     949              : END MODULE space_groups
        

Generated by: LCOV version 2.0-1