LCOV - code coverage report
Current view: top level - src/motion - space_groups.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 330 366 90.2 %
Date: 2024-04-20 06:29:22 Functions: 11 11 100.0 %

          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 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           8 :    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           8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: tmp_types
     102             :       LOGICAL                                            :: spglib
     103           8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: tmp_coor
     104             :       TYPE(spgr_type), POINTER                           :: spgr
     105             : 
     106           8 :       CALL timeset(routineN, handle)
     107             : 
     108           8 :       spgr => gopt_env%spgr
     109           8 :       CPASSERT(ASSOCIATED(spgr))
     110             : 
     111             :       !..total number of particles (atoms plus shells)
     112           8 :       IF (PRESENT(nparticle)) THEN
     113           8 :          CPASSERT(nparticle == SIZE(scoor, 2))
     114           8 :          spgr%nparticle = nparticle
     115             :       ELSE
     116           0 :          spgr%nparticle = SIZE(scoor, 2)
     117             :       END IF
     118             : 
     119           8 :       IF (PRESENT(n_atom)) THEN
     120           8 :          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           8 :       IF (PRESENT(n_core)) THEN
     130           8 :          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           8 :       IF (PRESENT(n_shell)) THEN
     136           8 :          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           8 :       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           8 :       spgr%nparticle_sym = spgr%nparticle
     146           8 :       spgr%n_atom_sym = spgr%n_atom
     147           8 :       spgr%n_core_sym = spgr%n_core
     148           8 :       spgr%n_shell_sym = spgr%n_shell
     149             : 
     150           8 :       spgr%iunit = iunit
     151           8 :       spgr%print_atoms = print_atoms
     152             : 
     153             :       ! accuracy for symmetry
     154           8 :       IF (PRESENT(eps_symmetry)) THEN
     155           8 :          spgr%eps_symmetry = eps_symmetry
     156             :       END IF
     157             : 
     158             :       ! vector to test reduced symmetry
     159           8 :       IF (PRESENT(pol)) THEN
     160           8 :          spgr%pol(1) = pol(1)
     161           8 :          spgr%pol(2) = pol(2)
     162           8 :          spgr%pol(3) = pol(3)
     163             :       END IF
     164             : 
     165          24 :       ALLOCATE (spgr%lat(spgr%nparticle))
     166         130 :       spgr%lat = .TRUE.
     167             : 
     168           8 :       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          40 :       ALLOCATE (tmp_coor(3, spgr%n_atom_sym), tmp_types(spgr%n_atom_sym))
     187             : 
     188           8 :       j = 0
     189          82 :       DO i = 1, spgr%n_atom
     190          82 :          IF (spgr%lat(i)) THEN
     191          74 :             j = j + 1
     192         296 :             tmp_coor(:, j) = scoor(:, i)
     193          74 :             tmp_types(j) = types(i)
     194             :          END IF
     195             :       END DO
     196             : 
     197             :       !..set cell values
     198           8 :       NULLIFY (spgr%cell_ref)
     199           8 :       CALL cell_create(spgr%cell_ref)
     200           8 :       CALL cell_copy(cell, spgr%cell_ref, tag="CELL_OPT_REF")
     201           8 :       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          16 :          SELECT CASE (gopt_env%cell_method_id)
     206             :          CASE (default_cell_direct_id)
     207           8 :             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           8 :             CPABORT("SPACE_GROUP_SYMMETRY invoked with an unknown optimization method.")
     214             :          END SELECT
     215             :       CASE DEFAULT
     216           8 :          CPABORT("SPACE_GROUP_SYMMETRY is not compatible with md.")
     217             :       END SELECT
     218             : 
     219             :       ! atom types
     220          24 :       ALLOCATE (spgr%atype(spgr%nparticle))
     221         130 :       spgr%atype(1:spgr%nparticle) = types(1:spgr%nparticle)
     222             : 
     223           8 :       spgr%n_operations = 0
     224             : 
     225             : #ifdef __SPGLIB
     226           8 :       spglib = .TRUE.
     227           8 :       CALL cite_reference(Togo2018)
     228             :       spgr%space_group_number = spg_get_international(spgr%international_symbol, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
     229         104 :                                                       spgr%n_atom_sym, eps_symmetry)
     230           8 :       buffer = ''
     231           8 :       nchars = strlcpy_c2f(buffer, spgr%international_symbol)
     232           8 :       spgr%international_symbol = buffer(1:nchars)
     233           8 :       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           8 :                                     spgr%n_atom_sym, eps_symmetry)
     239          40 :          ALLOCATE (spgr%rotations(3, 3, nop), spgr%translations(3, nop))
     240          32 :          ALLOCATE (spgr%eqatom(nop, spgr%nparticle))
     241          24 :          ALLOCATE (spgr%lop(nop))
     242           8 :          spgr%n_operations = nop
     243         420 :          spgr%lop = .TRUE.
     244             :          ierr = spg_get_symmetry(spgr%rotations, spgr%translations, nop, TRANSPOSE(cell%hmat), &
     245         104 :                                  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         104 :                                     spgr%n_atom_sym, eps_symmetry)
     249           8 :          buffer = ''
     250           8 :          nchars = strlcpy_c2f(buffer, spgr%schoenflies)
     251           8 :          spgr%schoenflies = buffer(1:nchars)
     252             : 
     253             :          ! Point Group
     254           8 :          tra_mat = 0
     255             :          ierr = spg_get_pointgroup(spgr%pointgroup_symbol, tra_mat, &
     256           8 :                                    spgr%rotations, spgr%n_operations)
     257           8 :          buffer = ''
     258           8 :          nchars = strlcpy_c2f(buffer, spgr%pointgroup_symbol)
     259           8 :          spgr%pointgroup_symbol = buffer(1:nchars)
     260             :       END IF
     261             : #else
     262             :       CPABORT("Symmetry library SPGLIB not available")
     263             :       spglib = .FALSE.
     264             : #endif
     265           8 :       spgr%symlib = spglib
     266             : 
     267           8 :       DEALLOCATE (tmp_coor, tmp_types)
     268             : 
     269           8 :       CALL timestop(handle)
     270             : 
     271           8 :    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           8 :    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           8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atype
     297           8 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: ranges
     298           8 :       INTEGER, DIMENSION(:), POINTER                     :: tmp
     299             :       LOGICAL                                            :: print_atoms
     300             :       REAL(KIND=dp)                                      :: eps_symmetry
     301           8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: scoord
     302           8 :       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           8 :       CALL timeset(routineN, handle)
     309             : 
     310           8 :       n_sr_rep = 0
     311           8 :       nparticle = 0
     312             :       n_atom = 0
     313           8 :       n_core = 0
     314           8 :       n_shell = 0
     315             : 
     316           8 :       NULLIFY (particles)
     317           8 :       NULLIFY (core_particles)
     318           8 :       NULLIFY (shell_particles)
     319             : 
     320             :       NULLIFY (cell)
     321           8 :       cell => subsys%cell
     322           8 :       CPASSERT(ASSOCIATED(cell))
     323             : 
     324           8 :       CALL cp_subsys_get(subsys, particles=particles, shell_particles=shell_particles, core_particles=core_particles)
     325             : 
     326           8 :       CPASSERT(ASSOCIATED(particles))
     327           8 :       n_atom = particles%n_els
     328             :       ! Check if we have other kinds of particles in this subsystem
     329           8 :       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           4 :       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           8 :       nparticle = n_atom + n_shell
     344          40 :       ALLOCATE (scoord(3, nparticle), atype(nparticle))
     345          82 :       DO i = 1, n_atom
     346          74 :          shell_index = particles%els(i)%shell_index
     347          82 :          IF (shell_index == 0) THEN
     348          26 :             CALL real_to_scaled(scoord(1:3, i), particles%els(i)%r(1:3), cell)
     349          26 :             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           8 :       CALL section_vals_val_get(geo_section, "SPGR_PRINT_ATOMS", l_val=print_atoms)
     360           8 :       CALL section_vals_val_get(geo_section, "EPS_SYMMETRY", r_val=eps_symmetry)
     361           8 :       CALL section_vals_val_get(geo_section, "SYMM_REDUCTION", r_vals=pol)
     362           8 :       CALL section_vals_val_get(geo_section, "SYMM_EXCLUDE_RANGE", n_rep_val=n_sr_rep)
     363           8 :       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           8 :                           n_core=n_core, n_shell=n_shell, iunit=iunit, print_atoms=print_atoms)
     377             :       END IF
     378             : 
     379             :       NULLIFY (spgr)
     380           8 :       spgr => gopt_env%spgr
     381             : 
     382           8 :       CALL spgr_find_equivalent_atoms(spgr, scoord)
     383           8 :       CALL spgr_reduce_symm(spgr)
     384           8 :       CALL spgr_rotations_subset(spgr)
     385             : 
     386           8 :       DEALLOCATE (scoord, atype)
     387             : 
     388           8 :       CALL timestop(handle)
     389             : 
     390          32 :    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           8 :    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           8 :       CALL timeset(routineN, handle)
     415             : 
     416           8 :       nop = spgr%n_operations
     417           8 :       natom = spgr%n_atom
     418           8 :       nshell = spgr%n_shell
     419             : 
     420           8 :       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         130 :       DO ia = 1, spgr%nparticle
     425        2158 :          spgr%eqatom(:, ia) = ia
     426             :       END DO
     427             : 
     428           8 :       !$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           8 :       !$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           8 :       CALL timestop(handle)
     483             : 
     484           8 :    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           8 :    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           8 :       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           8 :       CALL timeset(routineN, handle)
     506             : 
     507           8 :       nop = spgr%n_operations
     508           8 :       nparticle = spgr%nparticle
     509          40 :       ALLOCATE (x(3*nparticle), xold(3*nparticle))
     510         374 :       x = 0.0_dp
     511         130 :       DO ia = 1, nparticle
     512         122 :          ja = 3*(ia - 1)
     513         122 :          x(ja + 1) = x(ja + 1) + spgr%pol(1)
     514         122 :          x(ja + 2) = x(ja + 2) + spgr%pol(2)
     515         130 :          x(ja + 3) = x(ja + 3) + spgr%pol(3)
     516             :       END DO
     517         374 :       xold(:) = x(:)
     518             : 
     519             :       nops = 0
     520         420 :       DO ir = 1, nop
     521        6496 :          x = 0.d0
     522         412 :          spgr%lop(ir) = .TRUE.
     523        5356 :          rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
     524        2440 :          DO ia = 1, nparticle
     525        2028 :             IF (.NOT. spgr%lat(ia)) CYCLE
     526        2028 :             ja = 3*(ia - 1)
     527        8112 :             ri(1:3) = xold(ja + 1:ja + 3)
     528        2028 :             ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3)
     529        2028 :             ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3)
     530        2028 :             ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3)
     531        8524 :             x(ja + 1:ja + 3) = ro(1:3)
     532             :          END DO
     533        2440 :          DO ia = 1, nparticle
     534        2028 :             IF (.NOT. spgr%lat(ia)) CYCLE
     535        2028 :             ib = spgr%eqatom(ir, ia)
     536        2028 :             ja = 3*(ia - 1)
     537        2028 :             jb = 3*(ib - 1)
     538        8112 :             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        2440 :                             .AND. (ABS(ro(3)) < spgr%eps_symmetry))
     542             :          END DO
     543         420 :          IF (spgr%lop(ir)) nops = nops + 1
     544             :       END DO
     545             : 
     546           8 :       spgr%n_reduced_operations = nops
     547             : 
     548           8 :       DEALLOCATE (x, xold)
     549           8 :       CALL timestop(handle)
     550             : 
     551           8 :    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           8 :    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           8 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: mask
     570             : 
     571           8 :       CALL timeset(routineN, handle)
     572             : 
     573          24 :       ALLOCATE (mask(spgr%n_operations))
     574         420 :       mask = .TRUE.
     575             : 
     576         420 :       DO i = 1, spgr%n_operations
     577         420 :          IF (.NOT. spgr%lop(i)) mask(i) = .FALSE.
     578             :       END DO
     579             : 
     580         412 :       DO i = 1, spgr%n_operations - 1
     581         404 :          IF (.NOT. mask(i)) CYCLE
     582       16260 :          DO j = i + 1, spgr%n_operations
     583       16134 :             IF (.NOT. mask(j)) CYCLE
     584      121758 :             d(:, :) = spgr%rotations(:, :, j) - spgr%rotations(:, :, i)
     585      122162 :             IF (SUM(ABS(d)) == 0) mask(j) = .FALSE.
     586             :          END DO
     587             :       END DO
     588             : 
     589           8 :       spgr%n_operations_subset = 0
     590         420 :       DO i = 1, spgr%n_operations
     591         420 :          IF (mask(i)) spgr%n_operations_subset = spgr%n_operations_subset + 1
     592             :       END DO
     593             : 
     594          24 :       ALLOCATE (spgr%rotations_subset(3, 3, spgr%n_operations_subset))
     595             : 
     596           8 :       j = 0
     597         420 :       DO i = 1, spgr%n_operations
     598         420 :          IF (mask(i)) THEN
     599         124 :             j = j + 1
     600        1612 :             spgr%rotations_subset(:, :, j) = spgr%rotations(:, :, i)
     601             :          END IF
     602             :       END DO
     603             : 
     604           8 :       DEALLOCATE (mask)
     605           8 :       CALL timestop(handle)
     606             : 
     607           8 :    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           8 :    SUBROUTINE print_spgr(spgr)
     830             : 
     831             :       TYPE(spgr_type), INTENT(IN), POINTER               :: spgr
     832             : 
     833             :       INTEGER                                            :: i, j
     834             : 
     835           8 :       IF (spgr%iunit > 0) THEN
     836           4 :          WRITE (spgr%iunit, '(/,T2,A,A)') "----------------------------------------", &
     837           8 :             "---------------------------------------"
     838           4 :          WRITE (spgr%iunit, "(T2,A,T25,A,T77,A)") "----", "SPACE GROUP SYMMETRY INFORMATION", "----"
     839           4 :          WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
     840           8 :             "---------------------------------------"
     841           4 :          IF (spgr%symlib) THEN
     842           4 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| SPACE GROUP NUMBER:", &
     843           8 :                spgr%space_group_number
     844           4 :             WRITE (spgr%iunit, '(T2,A,T70,A11)') "SPGR| INTERNATIONAL SYMBOL:", &
     845           8 :                TRIM(ADJUSTR(spgr%international_symbol))
     846           4 :             WRITE (spgr%iunit, '(T2,A,T75,A6)') "SPGR| POINT GROUP SYMBOL:", &
     847           8 :                TRIM(ADJUSTR(spgr%pointgroup_symbol))
     848           4 :             WRITE (spgr%iunit, '(T2,A,T74,A7)') "SPGR| SCHOENFLIES SYMBOL:", &
     849           8 :                TRIM(ADJUSTR(spgr%schoenflies))
     850           4 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF SYMMETRY OPERATIONS:", &
     851           8 :                spgr%n_operations
     852           4 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF UNIQUE ROTATIONS:", &
     853           8 :                spgr%n_operations_subset
     854           4 :             WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF REDUCED SYMMETRY OPERATIONS:", &
     855           8 :                spgr%n_reduced_operations
     856           4 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF PARTICLES AND SYMMETRIZED PARTICLES:", &
     857           8 :                spgr%nparticle, spgr%nparticle_sym
     858           4 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF ATOMS AND SYMMETRIZED ATOMS:", &
     859           8 :                spgr%n_atom, spgr%n_atom_sym
     860           4 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF CORES AND SYMMETRIZED CORES:", &
     861           8 :                spgr%n_core, spgr%n_core_sym
     862           4 :             WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF SHELLS AND SYMMETRIZED SHELLS:", &
     863           8 :                spgr%n_shell, spgr%n_shell_sym
     864           4 :             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           8 :    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 1.15