LCOV - code coverage report
Current view: top level - src - molsym.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.8 % 657 636
Test Date: 2025-07-25 12:55:17 Functions: 96.8 % 31 30

            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 Molecular symmetry routines
      10              : !> \par History
      11              : !>      2008 adapted from older routines by Matthias Krack
      12              : !> \author jgh
      13              : ! **************************************************************************************************
      14              : MODULE molsym
      15              : 
      16              :    USE kinds,                           ONLY: dp
      17              :    USE mathconstants,                   ONLY: pi
      18              :    USE mathlib,                         ONLY: angle,&
      19              :                                               build_rotmat,&
      20              :                                               jacobi,&
      21              :                                               reflect_vector,&
      22              :                                               rotate_vector,&
      23              :                                               unit_matrix,&
      24              :                                               vector_product
      25              :    USE physcon,                         ONLY: angstrom
      26              :    USE util,                            ONLY: sort
      27              : #include "./base/base_uses.f90"
      28              : 
      29              :    IMPLICIT NONE
      30              :    PRIVATE
      31              : 
      32              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      33              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molsym'
      34              : 
      35              :    INTEGER, PARAMETER :: maxcn = 20, &
      36              :                          maxsec = maxcn + 1, &
      37              :                          maxses = 2*maxcn + 1, &
      38              :                          maxsig = maxcn + 1, &
      39              :                          maxsn = 2*maxcn
      40              : 
      41              :    PUBLIC :: molsym_type
      42              :    PUBLIC :: release_molsym, molecular_symmetry, print_symmetry
      43              : 
      44              : ! **************************************************************************************************
      45              : !> \brief Container for information about molecular symmetry
      46              : !> \param ain               : Atomic identification numbers (symmetry code).
      47              : !> \param aw                : Atomic weights for the symmetry analysis.
      48              : !> \param eps_geo           : Accuracy requested for analysis
      49              : !> \param inptostd          : Transformation matrix for input to standard orientation.
      50              : !> \param point_group_symbol: Point group symbol.
      51              : !> \param rotmat            : Rotation matrix.
      52              : !> \param sec               : List of C axes
      53              : !>                            (sec(:,i,j) => x,y,z of the ith j-fold C axis).
      54              : !> \param ses               : List of S axes
      55              : !>                            (ses(:,i,j) => x,y,z of the ith j-fold S axis).
      56              : !> \param sig               : List of mirror planes
      57              : !>                            (sig(:,i) => x,y,z of the ith mirror plane).
      58              : !> \param center_of_mass    : Shift vector for the center of mass.
      59              : !> \param tenmat            : Molecular tensor of inertia.
      60              : !> \param tenval            : Eigenvalues of the molecular tensor of inertia.
      61              : !> \param tenvec            : Eigenvectors of the molecular tensor of inertia.
      62              : !> \param group_of          : Group of equivalent atom.
      63              : !> \param llequatom         : Lower limit of a group in symequ_list.
      64              : !> \param ncn               : Degree of the C axis with highest degree.
      65              : !> \param ndod              : Number of substituted dodecahedral angles.
      66              : !> \param nequatom          : Number of equivalent atoms.
      67              : !> \param ngroup            : Number of groups of equivalent atoms.
      68              : !> \param nico              : Number of substituted icosahedral angles.
      69              : !> \param nlin              : Number of substituted angles of 180 degrees.
      70              : !> \param nsec              : Number of C axes.
      71              : !> \param nses              : Number of S axes.
      72              : !> \param nsig              : Number of mirror planes.
      73              : !> \param nsn               : Degree of the S axis with highest degree.
      74              : !> \param ntet              : Number of substituted tetrahedral angles.
      75              : !> \param point_group_order : Group order.
      76              : !> \param symequ_list       : List of all atoms ordered in groups of equivalent atoms.
      77              : !> \param ulequatom         : Upper limit of a group in symequ_list.
      78              : !> \param cubic : .TRUE., if a cubic point group was found (T,Th,Td,O,Oh,I,Ih).
      79              : !> \param dgroup: .TRUE., if a point group of D symmetry was found (Dn,Dnh,Dnd)
      80              : !> \param igroup: .TRUE., if a point group of icosahedral symmetry was found (I,Ih).
      81              : !> \param invers: .TRUE., if the molecule has a center of inversion.
      82              : !> \param linear: .TRUE., if the molecule is linear.
      83              : !> \param maxis : .TRUE., if the molecule has a main axis.
      84              : !> \param ogroup: .TRUE., if a point group of octahedral symmetry was found (O,Oh)
      85              : !> \param sgroup: .TRUE., if a point group of S symmetry was found.
      86              : !> \param sigmad: .TRUE., if there is a sigma_d mirror plane.
      87              : !> \param sigmah: .TRUE., if there is a sigma_h mirror plane.
      88              : !> \param sigmav: .TRUE., if there is a sigma_v mirror plane.
      89              : !> \param tgroup: .TRUE., if a point group of tetrahedral symmetry was found (T,Th,Td).
      90              : !> \par History
      91              : !>      05.2008 created
      92              : !> \author jgh
      93              : ! **************************************************************************************************
      94              :    TYPE molsym_type
      95              :       CHARACTER(LEN=4)                            :: point_group_symbol = ""
      96              :       INTEGER                                    :: point_group_order = -1
      97              :       INTEGER                                    :: ncn = -1, ndod = -1, ngroup = -1, nico = -1, &
      98              :                                                     nlin = -1, nsig = -1, nsn = -1, ntet = -1
      99              :       LOGICAL                                    :: cubic = .FALSE., dgroup = .FALSE., igroup = .FALSE., &
     100              :                                                     invers = .FALSE., linear = .FALSE., maxis = .FALSE., &
     101              :                                                     ogroup = .FALSE., sgroup = .FALSE., sigmad = .FALSE., &
     102              :                                                     sigmah = .FALSE., sigmav = .FALSE., tgroup = .FALSE.
     103              :       REAL(KIND=dp)                              :: eps_geo = 0.0_dp
     104              :       REAL(KIND=dp), DIMENSION(3)                :: center_of_mass = 0.0_dp, tenval = 0.0_dp
     105              :       REAL(KIND=dp), DIMENSION(3)                :: x_axis = 0.0_dp, y_axis = 0.0_dp, z_axis = 0.0_dp
     106              :       REAL(KIND=dp), DIMENSION(:), POINTER       :: ain => NULL(), aw => NULL()
     107              :       REAL(KIND=dp), DIMENSION(3, 3)              :: inptostd = 0.0_dp, rotmat = 0.0_dp, tenmat = 0.0_dp, tenvec = 0.0_dp
     108              :       REAL(KIND=dp), DIMENSION(3, maxsig)         :: sig = 0.0_dp
     109              :       REAL(KIND=dp), DIMENSION(3, maxsec, 2:maxcn) :: sec = 0.0_dp
     110              :       REAL(KIND=dp), DIMENSION(3, maxses, 2:maxsn) :: ses = 0.0_dp
     111              :       INTEGER, DIMENSION(maxcn)                  :: nsec = -1
     112              :       INTEGER, DIMENSION(maxsn)                  :: nses = -1
     113              :       INTEGER, DIMENSION(:), POINTER             :: group_of => NULL(), llequatom => NULL(), nequatom => NULL(), &
     114              :                                                     symequ_list => NULL(), ulequatom => NULL()
     115              :    END TYPE molsym_type
     116              : 
     117              : ! **************************************************************************************************
     118              : 
     119              : CONTAINS
     120              : 
     121              : ! **************************************************************************************************
     122              : !> \brief  Create an object of molsym type
     123              : !> \param sym ...
     124              : !> \param natoms ...
     125              : !> \author jgh
     126              : ! **************************************************************************************************
     127           58 :    SUBROUTINE create_molsym(sym, natoms)
     128              :       TYPE(molsym_type), POINTER                         :: sym
     129              :       INTEGER, INTENT(IN)                                :: natoms
     130              : 
     131           58 :       IF (ASSOCIATED(sym)) CALL release_molsym(sym)
     132              : 
     133       479718 :       ALLOCATE (sym)
     134              : 
     135              :       ALLOCATE (sym%ain(natoms), sym%aw(natoms), sym%group_of(natoms), sym%llequatom(natoms), &
     136          870 :                 sym%nequatom(natoms), sym%symequ_list(natoms), sym%ulequatom(natoms))
     137              : 
     138           58 :    END SUBROUTINE create_molsym
     139              : 
     140              : ! **************************************************************************************************
     141              : !> \brief  release an object of molsym type
     142              : !> \param sym ...
     143              : !> \author jgh
     144              : ! **************************************************************************************************
     145           58 :    SUBROUTINE release_molsym(sym)
     146              :       TYPE(molsym_type), POINTER                         :: sym
     147              : 
     148           58 :       CPASSERT(ASSOCIATED(sym))
     149              : 
     150           58 :       IF (ASSOCIATED(sym%aw)) THEN
     151           58 :          DEALLOCATE (sym%aw)
     152              :       END IF
     153           58 :       IF (ASSOCIATED(sym%ain)) THEN
     154           58 :          DEALLOCATE (sym%ain)
     155              :       END IF
     156           58 :       IF (ASSOCIATED(sym%group_of)) THEN
     157           58 :          DEALLOCATE (sym%group_of)
     158              :       END IF
     159           58 :       IF (ASSOCIATED(sym%llequatom)) THEN
     160           58 :          DEALLOCATE (sym%llequatom)
     161              :       END IF
     162           58 :       IF (ASSOCIATED(sym%nequatom)) THEN
     163           58 :          DEALLOCATE (sym%nequatom)
     164              :       END IF
     165           58 :       IF (ASSOCIATED(sym%symequ_list)) THEN
     166           58 :          DEALLOCATE (sym%symequ_list)
     167              :       END IF
     168           58 :       IF (ASSOCIATED(sym%ulequatom)) THEN
     169           58 :          DEALLOCATE (sym%ulequatom)
     170              :       END IF
     171              : 
     172           58 :       DEALLOCATE (sym)
     173              : 
     174           58 :    END SUBROUTINE release_molsym
     175              : 
     176              : ! **************************************************************************************************
     177              : 
     178              : ! **************************************************************************************************
     179              : !> \brief Add a new C_n axis to the list sec, but first check, if the
     180              : !>         Cn axis is already in the list.
     181              : !> \param n ...
     182              : !> \param a ...
     183              : !> \param sym ...
     184              : !> \par History
     185              : !>        Creation (19.10.98, Matthias Krack)
     186              : !> \author jgh
     187              : ! **************************************************************************************************
     188         1419 :    SUBROUTINE addsec(n, a, sym)
     189              :       INTEGER, INTENT(IN)                                :: n
     190              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: a
     191              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     192              : 
     193              :       INTEGER                                            :: isec
     194              :       REAL(dp)                                           :: length_of_a, scapro
     195              :       REAL(dp), DIMENSION(3)                             :: d
     196              : 
     197         1419 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
     198         5676 :       d(:) = a(:)/length_of_a
     199              : 
     200              :       ! Check, if the current Cn axis is already in the list
     201         6384 :       DO isec = 1, sym%nsec(n)
     202         6088 :          scapro = sym%sec(1, isec, n)*d(1) + sym%sec(2, isec, n)*d(2) + sym%sec(3, isec, n)*d(3)
     203         6384 :          IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
     204              :       END DO
     205          296 :       sym%ncn = MAX(sym%ncn, n)
     206              : 
     207              :       ! Add the new Cn axis to the list sec
     208          296 :       CPASSERT(sym%nsec(n) < maxsec)
     209          296 :       sym%nsec(1) = sym%nsec(1) + 1
     210          296 :       sym%nsec(n) = sym%nsec(n) + 1
     211         1184 :       sym%sec(:, sym%nsec(n), n) = d(:)
     212              : 
     213              :    END SUBROUTINE addsec
     214              : 
     215              : ! **************************************************************************************************
     216              : !> \brief  Add a new Sn axis to the list ses, but first check, if the
     217              : !>         Sn axis is already in the list.
     218              : !> \param n ...
     219              : !> \param a ...
     220              : !> \param sym ...
     221              : !> \par History
     222              : !>        Creation (19.10.98, Matthias Krack)
     223              : !> \author jgh
     224              : ! **************************************************************************************************
     225          284 :    SUBROUTINE addses(n, a, sym)
     226              :       INTEGER, INTENT(IN)                                :: n
     227              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: a
     228              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     229              : 
     230              :       INTEGER                                            :: ises
     231              :       REAL(dp)                                           :: length_of_a, scapro
     232              :       REAL(dp), DIMENSION(3)                             :: d
     233              : 
     234          284 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
     235         1136 :       d(:) = a(:)/length_of_a
     236              : 
     237              :       ! Check, if the current Sn axis is already in the list
     238         1051 :       DO ises = 1, sym%nses(n)
     239          888 :          scapro = sym%ses(1, ises, n)*d(1) + sym%ses(2, ises, n)*d(2) + sym%ses(3, ises, n)*d(3)
     240         1051 :          IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
     241              :       END DO
     242          163 :       sym%nsn = MAX(sym%nsn, n)
     243              : 
     244              :       ! Add the new Sn axis to the list ses
     245          163 :       CPASSERT(sym%nses(n) < maxses)
     246          163 :       sym%nses(1) = sym%nses(1) + 1
     247          163 :       sym%nses(n) = sym%nses(n) + 1
     248          652 :       sym%ses(:, sym%nses(n), n) = d(:)
     249              : 
     250              :    END SUBROUTINE addses
     251              : 
     252              : ! **************************************************************************************************
     253              : !> \brief  Add a new mirror plane to the list sig, but first check, if the
     254              : !>         normal vector of the mirror plane is already in the list.
     255              : !> \param a ...
     256              : !> \param sym ...
     257              : !> \par History
     258              : !>        Creation (19.10.98, Matthias Krack)
     259              : !> \author jgh
     260              : ! **************************************************************************************************
     261          597 :    SUBROUTINE addsig(a, sym)
     262              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: a
     263              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     264              : 
     265              :       INTEGER                                            :: isig
     266              :       REAL(dp)                                           :: length_of_a, scapro
     267              :       REAL(dp), DIMENSION(3)                             :: d
     268              : 
     269          597 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
     270         2388 :       d(:) = a(:)/length_of_a
     271              : 
     272              :       ! Check, if the normal vector of the current mirror plane is already in the list
     273         2426 :       DO isig = 1, sym%nsig
     274         2273 :          scapro = sym%sig(1, isig)*d(1) + sym%sig(2, isig)*d(2) + sym%sig(3, isig)*d(3)
     275         2426 :          IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) RETURN
     276              :       END DO
     277              : 
     278              :       ! Add the normal vector of the new mirror plane to the list sig
     279          153 :       CPASSERT(sym%nsig < maxsig)
     280          153 :       sym%nsig = sym%nsig + 1
     281          612 :       sym%sig(:, sym%nsig) = d(:)
     282              : 
     283              :    END SUBROUTINE addsig
     284              : 
     285              : ! **************************************************************************************************
     286              : !> \brief  Symmetry handling for only one atom.
     287              : !> \param sym ...
     288              : !> \par History
     289              : !>        Creation (19.10.98, Matthias Krack)
     290              : !> \author jgh
     291              : ! **************************************************************************************************
     292            1 :    SUBROUTINE atomsym(sym)
     293              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     294              : 
     295              : ! Set point group symbol
     296              : 
     297            1 :       sym%point_group_symbol = "R(3)"
     298              : 
     299              :       ! Set variables like D*h
     300            1 :       sym%linear = .TRUE.
     301            1 :       sym%invers = .TRUE.
     302            1 :       sym%dgroup = .TRUE.
     303            1 :       sym%sigmah = .TRUE.
     304              : 
     305            1 :    END SUBROUTINE atomsym
     306              : 
     307              : ! **************************************************************************************************
     308              : !> \brief Search for point groups with AXial SYMmetry (Cn,Cnh,Cnv,Dn,Dnh,Dnd,S2n).
     309              : !> \param coord ...
     310              : !> \param sym ...
     311              : !> \par History
     312              : !>        Creation (19.10.98, Matthias Krack)
     313              : !> \author jgh
     314              : ! **************************************************************************************************
     315           45 :    SUBROUTINE axsym(coord, sym)
     316              :       REAL(Kind=dp), DIMENSION(:, :), INTENT(inout)      :: coord
     317              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     318              : 
     319              :       INTEGER                                            :: iatom, isig, jatom, m, n, natoms, nx, &
     320              :                                                             ny, nz
     321              :       REAL(dp)                                           :: phi
     322              :       REAL(dp), DIMENSION(3)                             :: a, b, d
     323              : 
     324              : ! Find the correct main axis and rotate main axis to z axis
     325              : 
     326           45 :       IF ((sym%ncn == 2) .AND. (sym%nses(2) == 3)) THEN
     327            2 :          IF (sym%nsn == 4) THEN
     328              :             ! Special case: D2d
     329            0 :             phi = angle(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
     330            0 :             d(:) = vector_product(sym%ses(:, 1, sym%nsn), sym%z_axis(:))
     331            0 :             CALL rotate_molecule(phi, d(:), sym, coord)
     332              :          ELSE
     333              :             ! Special cases: D2 and D2h
     334            2 :             phi = 0.5_dp*pi
     335            2 :             nx = naxis(sym%x_axis(:), coord, sym)
     336            2 :             ny = naxis(sym%y_axis(:), coord, sym)
     337            2 :             nz = naxis(sym%z_axis(:), coord, sym)
     338            2 :             IF ((nx > ny) .AND. (nx > nz)) THEN
     339            0 :                CALL rotate_molecule(-phi, sym%y_axis(:), sym, coord)
     340            2 :             ELSE IF ((ny > nz) .AND. (ny > nx)) THEN
     341            0 :                CALL rotate_molecule(phi, sym%x_axis(:), sym, coord)
     342              :             END IF
     343              :          END IF
     344              :       ELSE
     345           43 :          phi = angle(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
     346           43 :          d(:) = vector_product(sym%sec(:, 1, sym%ncn), sym%z_axis(:))
     347           43 :          CALL rotate_molecule(phi, d(:), sym, coord)
     348              :       END IF
     349              : 
     350              :       ! Search for C2 axes perpendicular to the main axis
     351              :       ! Loop over all pairs of atoms (except on the z axis)
     352           45 :       natoms = SIZE(coord, 2)
     353          459 :       DO iatom = 1, natoms
     354         1656 :          a(:) = coord(:, iatom)
     355          459 :          IF ((ABS(a(1)) > sym%eps_geo) .OR. (ABS(a(2)) > sym%eps_geo)) THEN
     356          388 :             a(3) = 0.0_dp
     357          388 :             IF (caxis(2, a(:), sym, coord)) CALL addsec(2, a(:), sym)
     358          388 :             d(:) = vector_product(a(:), sym%z_axis(:))
     359          388 :             IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
     360         2283 :             DO jatom = iatom + 1, natoms
     361         7580 :                b(:) = coord(:, jatom)
     362         2309 :                IF ((ABS(b(1)) > sym%eps_geo) .OR. (ABS(b(2)) > sym%eps_geo)) THEN
     363         1871 :                   b(3) = 0.0_dp
     364         1871 :                   phi = 0.5_dp*angle(a(:), b(:))
     365         1871 :                   d(:) = vector_product(a(:), b(:))
     366         1871 :                   b(:) = rotate_vector(a(:), phi, d(:))
     367         1871 :                   IF (caxis(2, b(:), sym, coord)) CALL addsec(2, b(:), sym)
     368         1871 :                   d(:) = vector_product(b(:), sym%z_axis)
     369         1895 :                   IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
     370              :                END IF
     371              :             END DO
     372              :          END IF
     373              :       END DO
     374              : 
     375              :       ! Check the xy plane for mirror plane
     376           45 :       IF (sigma(sym%z_axis(:), sym, coord)) THEN
     377           14 :          CALL addsig(sym%z_axis(:), sym)
     378           14 :          sym%sigmah = .TRUE.
     379              :       END IF
     380              : 
     381              :       ! Set logical variables for point group determination ***
     382           45 :       IF (sym%nsec(2) > 1) THEN
     383           21 :          sym%dgroup = .TRUE.
     384           21 :          IF ((.NOT. sym%sigmah) .AND. (sym%nsig > 0)) sym%sigmad = .TRUE.
     385              :       ELSE
     386           24 :          IF ((.NOT. sym%sigmah) .AND. (sym%nsig > 0)) sym%sigmav = .TRUE.
     387              :          ! Cnh groups with n>3 were wrong
     388              :          ! IF (sym%nsn > 3) sym%sgroup = .TRUE.
     389           24 :          IF ((.NOT. sym%sigmah) .AND. (sym%nsn > 3)) sym%sgroup = .TRUE.
     390              :       END IF
     391              : 
     392              :       ! Rotate to standard orientation
     393           45 :       n = 0
     394           45 :       m = 0
     395           45 :       IF ((sym%dgroup .AND. sym%sigmah) .OR. (sym%dgroup .AND. sym%sigmad) .OR. sym%sigmav) THEN
     396              :          ! Dnh, Dnd or Cnv
     397          133 :          DO isig = 1, sym%nsig
     398          133 :             IF (ABS(sym%sig(3, isig)) < sym%eps_geo) THEN
     399          105 :                n = nsigma(sym%sig(:, isig), sym, coord)
     400          105 :                IF (n > m) THEN
     401           21 :                   m = n
     402           84 :                   a(:) = sym%sig(:, isig)
     403              :                END IF
     404              :             END IF
     405              :          END DO
     406           21 :          IF (m > 0) THEN
     407              :             ! Check for special case: C2v with all atoms in a plane
     408           21 :             IF (sym%sigmav .AND. (m == natoms)) THEN
     409            1 :                phi = angle(a(:), sym%x_axis(:))
     410            1 :                d(:) = vector_product(a(:), sym%x_axis(:))
     411              :             ELSE
     412           20 :                phi = angle(a(:), sym%y_axis(:))
     413           20 :                d(:) = vector_product(a(:), sym%y_axis(:))
     414              :             END IF
     415           21 :             CALL rotate_molecule(phi, d(:), sym, coord)
     416              :          END IF
     417           24 :       ELSE IF (sym%sigmah) THEN
     418              :          ! Cnh
     419           81 :          DO iatom = 1, natoms
     420           81 :             IF (ABS(coord(3, iatom)) < sym%eps_geo) THEN
     421           70 :                n = naxis(coord(:, iatom), coord, sym)
     422           70 :                IF (n > m) THEN
     423           28 :                   m = n
     424           28 :                   a(:) = coord(:, iatom)
     425              :                END IF
     426              :             END IF
     427              :          END DO
     428            7 :          IF (m > 0) THEN
     429            7 :             phi = angle(a(:), sym%x_axis(:))
     430            7 :             d(:) = vector_product(a(:), sym%x_axis(:))
     431            7 :             CALL rotate_molecule(phi, d(:), sym, coord)
     432              :          END IF
     433              :          ! No action for Dn, Cn or S2n ***
     434              :       END IF
     435              : 
     436           45 :    END SUBROUTINE axsym
     437              : 
     438              : ! **************************************************************************************************
     439              : !> \brief Generate a symmetry list to identify equivalent atoms.
     440              : !> \param sym ...
     441              : !> \param coord ...
     442              : !> \par History
     443              : !>        Creation (19.10.98, Matthias Krack)
     444              : !> \author jgh
     445              : ! **************************************************************************************************
     446           58 :    SUBROUTINE build_symequ_list(sym, coord)
     447              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     448              :       REAL(Kind=dp), DIMENSION(:, :), INTENT(inout)      :: coord
     449              : 
     450              :       INTEGER                                            :: i, iatom, icn, iequatom, incr, isec, &
     451              :                                                             ises, isig, isn, jatom, jcn, jsn, &
     452              :                                                             natoms
     453              :       REAL(KIND=dp)                                      :: phi
     454              :       REAL(KIND=dp), DIMENSION(3)                        :: a
     455              : 
     456           58 :       natoms = SIZE(coord, 2)
     457              : 
     458           58 :       IF (sym%sigmah) THEN
     459              :          incr = 1
     460              :       ELSE
     461           41 :          incr = 2
     462              :       END IF
     463              : 
     464              :       ! Initialize the atom and the group counter
     465           58 :       iatom = 1
     466           58 :       sym%ngroup = 1
     467              : 
     468              :       loop: DO
     469              : 
     470              :          ! Loop over all mirror planes
     471          332 :          DO isig = 1, sym%nsig
     472          226 :             a(:) = reflect_vector(coord(:, iatom), sym%sig(:, isig))
     473          226 :             iequatom = equatom(iatom, a(:), sym, coord)
     474          332 :             IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
     475          111 :                sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
     476          111 :                sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
     477          111 :                sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
     478              :             END IF
     479              :          END DO
     480              : 
     481              :          ! Loop over all Cn axes
     482          445 :          DO icn = 2, sym%ncn
     483          829 :             DO isec = 1, sym%nsec(icn)
     484         1465 :                DO jcn = 1, icn - 1
     485         1126 :                   IF (newse(jcn, icn)) THEN
     486          644 :                      phi = 2.0_dp*pi*REAL(jcn, KIND=dp)/REAL(icn, KIND=dp)
     487          644 :                      a(:) = rotate_vector(coord(:, iatom), phi, sym%sec(:, isec, icn))
     488          644 :                      iequatom = equatom(iatom, a(:), sym, coord)
     489          644 :                      IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
     490          328 :                         sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
     491          328 :                         sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
     492          328 :                         sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
     493              :                      END IF
     494              :                   END IF
     495              :                END DO
     496              :             END DO
     497              :          END DO
     498              : 
     499              :          ! Loop over all Sn axes
     500          339 :          DO isn = 2, sym%nsn
     501          442 :             DO ises = 1, sym%nses(isn)
     502          643 :                DO jsn = 1, isn - 1, incr
     503          410 :                   IF (newse(jsn, isn)) THEN
     504          243 :                      phi = 2.0_dp*pi*REAL(jsn, KIND=dp)/REAL(isn, KIND=dp)
     505          243 :                      a(:) = rotate_vector(coord(:, iatom), phi, sym%ses(:, ises, isn))
     506          972 :                      a(:) = reflect_vector(a(:), sym%ses(:, ises, isn))
     507          243 :                      iequatom = equatom(iatom, a(:), sym, coord)
     508          243 :                      IF ((iequatom > 0) .AND. (.NOT. in_symequ_list(iequatom, sym))) THEN
     509           30 :                         sym%ulequatom(sym%ngroup) = sym%ulequatom(sym%ngroup) + 1
     510           30 :                         sym%symequ_list(sym%ulequatom(sym%ngroup)) = iequatom
     511           30 :                         sym%nequatom(sym%ngroup) = sym%nequatom(sym%ngroup) + 1
     512              :                      END IF
     513              :                   END IF
     514              :                END DO
     515              :             END DO
     516              :          END DO
     517              : 
     518              :          ! Exit loop, if all atoms are in the list
     519          106 :          IF (sym%ulequatom(sym%ngroup) == natoms) EXIT
     520              : 
     521              :          ! Search for the next atom which is not in the list yet
     522          156 :          DO jatom = 2, natoms
     523           98 :             IF (.NOT. in_symequ_list(jatom, sym)) THEN
     524           48 :                iatom = jatom
     525           48 :                sym%ngroup = sym%ngroup + 1
     526           48 :                sym%llequatom(sym%ngroup) = sym%ulequatom(sym%ngroup - 1) + 1
     527           48 :                sym%ulequatom(sym%ngroup) = sym%llequatom(sym%ngroup)
     528           48 :                sym%symequ_list(sym%llequatom(sym%ngroup)) = iatom
     529           48 :                CYCLE loop
     530              :             END IF
     531              :          END DO
     532              : 
     533              :       END DO loop
     534              : 
     535              :       ! Generate list vector group_of
     536          164 :       DO i = 1, sym%ngroup
     537          739 :          DO iequatom = sym%llequatom(i), sym%ulequatom(i)
     538          681 :             sym%group_of(sym%symequ_list(iequatom)) = i
     539              :          END DO
     540              :       END DO
     541              : 
     542           58 :    END SUBROUTINE build_symequ_list
     543              : 
     544              : ! **************************************************************************************************
     545              : !> \brief Rotate the molecule about a n-fold axis defined by the vector a
     546              : !>        using a rotation angle of 2*pi/n. Check, if the axis is a Cn axis.
     547              : !> \param n ...
     548              : !> \param a ...
     549              : !> \param sym ...
     550              : !> \param coord ...
     551              : !> \return ...
     552              : !> \par History
     553              : !>        Creation (19.10.98, Matthias Krack)
     554              : !> \author jgh
     555              : ! **************************************************************************************************
     556         5543 :    FUNCTION caxis(n, a, sym, coord)
     557              :       INTEGER, INTENT(IN)                                :: n
     558              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
     559              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     560              :       REAL(Kind=dp), DIMENSION(:, :), INTENT(inout)      :: coord
     561              :       LOGICAL                                            :: caxis
     562              : 
     563              :       INTEGER                                            :: iatom, natoms
     564              :       REAL(KIND=dp)                                      :: length_of_a, phi
     565              :       REAL(KIND=dp), DIMENSION(3)                        :: b
     566              : 
     567         5543 :       caxis = .FALSE.
     568         5543 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
     569              : 
     570              :       ! Check the length of the axis vector a
     571         5543 :       natoms = SIZE(coord, 2)
     572         5543 :       IF (length_of_a > sym%eps_geo) THEN
     573              :          ! Calculate the rotation angle phi and build up the rotation matrix rotmat
     574         5409 :          phi = 2.0_dp*pi/REAL(n, dp)
     575         5409 :          CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
     576              :          ! Check all atoms
     577        18830 :          DO iatom = 1, natoms
     578              :             ! Rotate the actual atom by 2*pi/n about a
     579       234507 :             b(:) = MATMUL(sym%rotmat(:, :), coord(:, iatom))
     580              :             ! Check if the coordinates of the rotated atom
     581              :             ! are in the coordinate set of the molecule
     582        18830 :             IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
     583              :          END DO
     584              :          ! The axis defined by vector a is a Cn axis, thus return with caxis = .TRUE.
     585              :          caxis = .TRUE.
     586              :       END IF
     587              : 
     588              :    END FUNCTION caxis
     589              : 
     590              : ! **************************************************************************************************
     591              : !> \brief Check for a point group of cubic symmetry (I,Ih,O,Oh,T,Th,Td).
     592              : !> \param sym ...
     593              : !> \param coord ...
     594              : !> \param failed ...
     595              : !> \par History
     596              : !>        Creation (19.10.98, Matthias Krack)
     597              : !> \author jgh
     598              : ! **************************************************************************************************
     599            7 :    SUBROUTINE cubsym(sym, coord, failed)
     600              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     601              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
     602              :       LOGICAL, INTENT(INOUT)                             :: failed
     603              : 
     604              :       INTEGER                                            :: i, iatom, iax, ic3, isec, jatom, jc3, &
     605              :                                                             jsec, katom, natoms, nc3
     606              :       REAL(KIND=dp)                                      :: phi1, phi2, phidd, phidi, phiii
     607              :       REAL(KIND=dp), DIMENSION(3)                        :: a, b, c, d
     608              : 
     609              :       ! Angle between two adjacent atoms of the dodecahedron => <(C3,C3)
     610            7 :       phidd = ATAN(0.4_dp*SQRT(5.0_dp))
     611              :       ! Angle between two adjacent atoms of the icosahedron and the dodecahedron => <(C5,C3)
     612            7 :       phidi = ATAN(3.0_dp - SQRT(5.0_dp))
     613              :       ! Angle between two adjacent atoms of the icosahedron <(C5,C5)
     614            7 :       phiii = ATAN(2.0_dp)
     615              : 
     616              :       ! Find a pair of C3 axes
     617            7 :       natoms = SIZE(coord, 2)
     618          146 :       DO iatom = 1, natoms
     619              :          ! Check all atomic vectors for C3 axis
     620          146 :          IF (caxis(3, coord(:, iatom), sym, coord)) THEN
     621            2 :             CALL addsec(3, coord(:, iatom), sym)
     622            2 :             IF (sym%nsec(3) > 1) EXIT
     623              :          END IF
     624              :       END DO
     625              : 
     626              :       ! Check the center of mass vector of a triangle
     627              :       ! defined by three atomic vectors for C3 axis
     628            7 :       IF (sym%nsec(3) < 2) THEN
     629              : 
     630            6 :          loop: DO iatom = 1, natoms
     631           24 :             a(:) = coord(:, iatom)
     632           24 :             DO jatom = iatom + 1, natoms
     633          734 :                DO katom = jatom + 1, natoms
     634              :                   IF ((ABS(dist(coord(:, iatom), coord(:, jatom)) &
     635          716 :                            - dist(coord(:, iatom), coord(:, katom))) < sym%eps_geo) .AND. &
     636              :                       (ABS(dist(coord(:, iatom), coord(:, jatom)) &
     637           18 :                            - dist(coord(:, jatom), coord(:, katom))) < sym%eps_geo)) THEN
     638           48 :                      b(:) = a(:) + coord(:, jatom) + coord(:, katom)
     639           12 :                      IF (caxis(3, b(:), sym, coord)) THEN
     640           12 :                         CALL addsec(3, b(:), sym)
     641           12 :                         IF (sym%nsec(3) > 1) EXIT loop
     642              :                      END IF
     643              :                   END IF
     644              :                END DO
     645              :             END DO
     646              :          END DO loop
     647              : 
     648              :       END IF
     649              : 
     650              :       ! Return after unsuccessful search for a pair of C3 axes
     651            7 :       IF (sym%nsec(3) < 2) THEN
     652            0 :          failed = .TRUE.
     653            0 :          RETURN
     654              :       END IF
     655              : 
     656              :       ! Generate the remaining C3 axes
     657           16 :       nc3 = 0
     658              :       DO
     659           16 :          nc3 = sym%nsec(3)
     660           82 :          DO ic3 = 1, nc3
     661          264 :             a(:) = sym%sec(:, ic3, 3)
     662          462 :             DO jc3 = 1, nc3
     663          446 :                IF (ic3 /= jc3) THEN
     664         1256 :                   d(:) = sym%sec(:, jc3, 3)
     665          942 :                   DO i = 1, 2
     666          628 :                      phi1 = 2.0_dp*REAL(i, KIND=dp)*pi/3.0_dp
     667          628 :                      b(:) = rotate_vector(a(:), phi1, d(:))
     668          942 :                      CALL addsec(3, b(:), sym)
     669              :                   END DO
     670              :                END IF
     671              :             END DO
     672              :          END DO
     673              : 
     674           16 :          IF (sym%nsec(3) > nc3) THEN
     675              :             ! New C3 axes were found
     676              :             CYCLE
     677              :          ELSE
     678              :             ! In the case of I or Ih there have to be more C3 axes
     679            7 :             IF (sym%nsec(3) == 4) THEN
     680           20 :                a(:) = sym%sec(:, 1, 3)
     681           20 :                b(:) = sym%sec(:, 2, 3)
     682            5 :                phi1 = 0.5_dp*angle(a(:), b(:))
     683            5 :                IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
     684            5 :                d(:) = vector_product(a(:), b(:))
     685            5 :                b(:) = rotate_vector(a(:), phi1, d(:))
     686           20 :                c(:) = sym%sec(:, 3, 3)
     687            5 :                phi1 = 0.5_dp*angle(a(:), c(:))
     688            5 :                IF (phi1 < 0.5_dp*pi) phi1 = phi1 - 0.5_dp*pi
     689            5 :                d(:) = vector_product(a(:), c(:))
     690            5 :                c(:) = rotate_vector(a(:), phi1, d(:))
     691            5 :                d(:) = vector_product(b(:), c(:))
     692            5 :                phi1 = 0.5_dp*phidd
     693            5 :                a(:) = rotate_vector(b(:), phi1, d(:))
     694            5 :                IF (caxis(3, a(:), sym, coord)) THEN
     695            0 :                   CALL addsec(3, a(:), sym)
     696              :                ELSE
     697            5 :                   phi2 = 0.5_dp*pi - phi1
     698            5 :                   a(:) = rotate_vector(b(:), phi2, d(:))
     699            5 :                   IF (caxis(3, a(:), sym, coord)) CALL addsec(3, a(:), sym)
     700              :                END IF
     701              : 
     702            5 :                IF (sym%nsec(3) > 4) CYCLE
     703              : 
     704            2 :             ELSE IF (sym%nsec(3) /= 10) THEN
     705              : 
     706              :                !  C3 axes found, but only 4 or 10 are possible
     707            0 :                failed = .TRUE.
     708            0 :                RETURN
     709              : 
     710              :             END IF
     711              : 
     712              :             ! Exit loop, if 4 or 10 C3 axes were found
     713            9 :             EXIT
     714              : 
     715              :          END IF
     716              : 
     717              :       END DO
     718              : 
     719              :       ! Loop over all pairs of C3 axes to find all C5, C4 and C2 axes
     720           47 :       DO isec = 1, sym%nsec(3)
     721              : 
     722          160 :          a(:) = sym%sec(:, isec, 3)
     723          167 :          DO jsec = isec + 1, sym%nsec(3)
     724          120 :             phi1 = 0.5_dp*angle(a(:), sym%sec(:, jsec, 3))
     725          120 :             d(:) = vector_product(a(:), sym%sec(:, jsec, 3))
     726              : 
     727              :             ! Check for C5 axis (I,Ih)
     728          120 :             IF (sym%nsec(3) == 10) THEN
     729           90 :                b(:) = rotate_vector(a(:), phidi, d(:))
     730           90 :                IF (caxis(5, b(:), sym, coord)) THEN
     731           14 :                   CALL addsec(5, b(:), sym)
     732           14 :                   phi1 = phidi + phiii
     733           14 :                   b(:) = rotate_vector(a(:), phi1, d(:))
     734           14 :                   IF (caxis(5, b(:), sym, coord)) CALL addsec(5, b(:), sym)
     735              :                END IF
     736              :             END IF
     737              : 
     738              :             ! Check for C4 (O,Oh), C2 and S2 (center of inversion) axis
     739          360 :             DO i = 0, 1
     740          240 :                phi2 = phi1 - 0.5_dp*REAL(i, KIND=dp)*pi
     741          240 :                b(:) = rotate_vector(a(:), phi2, d(:))
     742          240 :                IF (caxis(2, b(:), sym, coord)) THEN
     743          134 :                   CALL addsec(2, b(:), sym)
     744          134 :                   IF (sym%nsec(3) == 4) THEN
     745           42 :                      IF (caxis(4, b(:), sym, coord)) CALL addsec(4, b(:), sym)
     746              :                   END IF
     747          134 :                   IF (saxis(2, b(:), sym, coord)) THEN
     748           64 :                      CALL addses(2, b(:), sym)
     749           64 :                      sym%invers = .TRUE.
     750              :                   END IF
     751              :                END IF
     752          360 :                IF (sigma(b(:), sym, coord)) CALL addsig(b(:), sym)
     753              :             END DO
     754              : 
     755              :             ! Check for mirror plane
     756          160 :             IF (sigma(d(:), sym, coord)) CALL addsig(d(:), sym)
     757              : 
     758              :          END DO
     759              : 
     760              :       END DO
     761              : 
     762              :       ! Set the logical variables for point group determination
     763            7 :       iax = sym%ncn
     764              : 
     765            7 :       IF ((sym%ncn == 5) .AND. (sym%nsec(5) > 1)) THEN
     766            2 :          sym%igroup = .TRUE.
     767            5 :       ELSE IF ((sym%ncn == 4) .AND. (sym%nsec(4) > 1)) THEN
     768            2 :          sym%ogroup = .TRUE.
     769            3 :       ELSE IF ((sym%ncn == 3) .AND. (sym%nsec(3) > 1)) THEN
     770            3 :          sym%tgroup = .TRUE.
     771            3 :          IF ((.NOT. sym%invers) .AND. (sym%nsig > 0)) sym%sigmad = .TRUE.
     772              :          iax = 2
     773              :       ELSE
     774              :          ! No main axis found
     775            0 :          failed = .TRUE.
     776            0 :          RETURN
     777              :       END IF
     778              : 
     779              :       ! Rotate molecule to standard orientation
     780            7 :       phi1 = angle(sym%sec(:, 1, iax), sym%z_axis(:))
     781            7 :       d(:) = vector_product(sym%sec(:, 1, iax), sym%z_axis(:))
     782            7 :       CALL rotate_molecule(phi1, d(:), sym, coord)
     783              : 
     784           28 :       a(:) = sym%sec(:, 2, iax)
     785            7 :       a(3) = 0.0_dp
     786            7 :       phi2 = angle(a(:), sym%x_axis(:))
     787            7 :       d(:) = vector_product(a(:), sym%x_axis(:))
     788            7 :       CALL rotate_molecule(phi2, d(:), sym, coord)
     789              : 
     790              :    END SUBROUTINE cubsym
     791              : 
     792              : ! **************************************************************************************************
     793              : !> \brief The number of a equivalent atom to atoma is returned. If there
     794              : !>        is no equivalent atom, then zero is returned. The cartesian
     795              : !>        coordinates of the equivalent atom are defined by the vector a.
     796              : !> \param atoma ...
     797              : !> \param a ...
     798              : !> \param sym ...
     799              : !> \param coord ...
     800              : !> \return ...
     801              : !> \par History
     802              : !>        Creation (19.10.98, Matthias Krack)
     803              : !> \author jgh
     804              : ! **************************************************************************************************
     805        38170 :    FUNCTION equatom(atoma, a, sym, coord)
     806              :       INTEGER, INTENT(IN)                                :: atoma
     807              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
     808              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     809              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
     810              :       INTEGER                                            :: equatom
     811              : 
     812              :       INTEGER                                            :: iatom, natoms
     813              : 
     814        38170 :       equatom = 0
     815        38170 :       natoms = SIZE(coord, 2)
     816       413408 :       DO iatom = 1, natoms
     817              :          IF ((ABS(sym%ain(iatom) - sym%ain(atoma)) < TINY(0.0_dp)) .AND. &
     818              :              (ABS(a(1) - coord(1, iatom)) < sym%eps_geo) .AND. &
     819       400638 :              (ABS(a(2) - coord(2, iatom)) < sym%eps_geo) .AND. &
     820        12770 :              (ABS(a(3) - coord(3, iatom)) < sym%eps_geo)) THEN
     821        38170 :             equatom = iatom
     822              :             RETURN
     823              :          END IF
     824              :       END DO
     825              : 
     826              :    END FUNCTION equatom
     827              : 
     828              : ! **************************************************************************************************
     829              : !> \brief Calculate the order of the point group.
     830              : !> \param sym ...
     831              : !> \par History
     832              : !>        Creation (19.10.98, Matthias Krack)
     833              : !> \author jgh
     834              : ! **************************************************************************************************
     835           55 :    SUBROUTINE get_point_group_order(sym)
     836              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     837              : 
     838              :       INTEGER                                            :: icn, incr, isec, ises, isn, jcn, jsn
     839              : 
     840              :       ! Count all symmetry elements of the symmetry group
     841              :       ! First E and all mirror planes
     842           55 :       sym%point_group_order = 1 + sym%nsig
     843              : 
     844              :       ! Loop over all C axes
     845          249 :       DO icn = 2, sym%ncn
     846          545 :          DO isec = 1, sym%nsec(icn)
     847         1021 :             DO jcn = 1, icn - 1
     848          827 :                IF (newse(jcn, icn)) sym%point_group_order = sym%point_group_order + 1
     849              :             END DO
     850              :          END DO
     851              :       END DO
     852              : 
     853              :       ! Loop over all S axes
     854           55 :       IF (sym%sigmah) THEN
     855              :          incr = 1
     856              :       ELSE
     857           40 :          incr = 2
     858              :       END IF
     859              : 
     860          212 :       DO isn = 2, sym%nsn
     861          285 :          DO ises = 1, sym%nses(isn)
     862          452 :             DO jsn = 1, isn - 1, incr
     863          295 :                IF (newse(jsn, isn)) sym%point_group_order = sym%point_group_order + 1
     864              :             END DO
     865              :          END DO
     866              :       END DO
     867              : 
     868           55 :    END SUBROUTINE get_point_group_order
     869              : 
     870              : ! **************************************************************************************************
     871              : !> \brief Get the point group symbol.
     872              : !> \param sym ...
     873              : !> \par History
     874              : !>        Creation (19.10.98, Matthias Krack)
     875              : !> \author jgh
     876              : ! **************************************************************************************************
     877           57 :    SUBROUTINE get_point_group_symbol(sym)
     878              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     879              : 
     880              :       CHARACTER(LEN=4)                                   :: degree
     881              : 
     882           57 :       IF (sym%linear) THEN
     883            2 :          IF (sym%invers) THEN
     884            1 :             sym%point_group_symbol = "D*h"
     885              :          ELSE
     886            1 :             sym%point_group_symbol = "C*v"
     887              :          END IF
     888           55 :       ELSE IF (sym%cubic) THEN
     889            7 :          IF (sym%igroup) THEN
     890            2 :             sym%point_group_symbol = "I"
     891            5 :          ELSE IF (sym%ogroup) THEN
     892            2 :             sym%point_group_symbol = "O"
     893            3 :          ELSE IF (sym%tgroup) THEN
     894            3 :             sym%point_group_symbol = "T"
     895              :          END IF
     896            7 :          IF (sym%invers) THEN
     897            3 :             sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
     898            4 :          ELSE IF (sym%sigmad) THEN
     899            1 :             sym%point_group_symbol = TRIM(sym%point_group_symbol)//"d"
     900              :          END IF
     901           48 :       ELSE IF (sym%maxis) THEN
     902           45 :          IF (sym%dgroup) THEN
     903           21 :             WRITE (degree, "(I3)") sym%ncn
     904           21 :             sym%point_group_symbol = "D"//TRIM(ADJUSTL(degree))
     905           21 :             IF (sym%sigmah) THEN
     906            7 :                sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
     907           14 :             ELSE IF (sym%sigmad) THEN
     908            7 :                sym%point_group_symbol = TRIM(sym%point_group_symbol)//"d"
     909              :             END IF
     910           24 :          ELSE IF (sym%sgroup) THEN
     911            3 :             WRITE (degree, "(I3)") sym%nsn
     912            3 :             sym%point_group_symbol = "S"//TRIM(ADJUSTL(degree))
     913              :          ELSE
     914           21 :             WRITE (degree, "(I3)") sym%ncn
     915           21 :             sym%point_group_symbol = "C"//TRIM(ADJUSTL(degree))
     916           21 :             IF (sym%sigmah) THEN
     917            7 :                sym%point_group_symbol = TRIM(sym%point_group_symbol)//"h"
     918           14 :             ELSE IF (sym%sigmav) THEN
     919            7 :                sym%point_group_symbol = TRIM(sym%point_group_symbol)//"v"
     920              :             END IF
     921              :          END IF
     922              :       ELSE
     923            3 :          IF (sym%invers) THEN
     924            1 :             sym%point_group_symbol = "Ci"
     925            2 :          ELSE IF (sym%sigmah) THEN
     926            1 :             sym%point_group_symbol = "Cs"
     927              :          ELSE
     928            1 :             sym%point_group_symbol = "C1"
     929              :          END IF
     930              :       END IF
     931              : 
     932           57 :    END SUBROUTINE get_point_group_symbol
     933              : 
     934              : ! **************************************************************************************************
     935              : !> \brief Initialization of the global variables of module symmetry.
     936              : !> \param sym ...
     937              : !> \param atype ...
     938              : !> \param weight ...
     939              : !> \par History
     940              : !>        Creation (19.10.98, Matthias Krack)
     941              : !> \author jgh
     942              : ! **************************************************************************************************
     943           58 :    SUBROUTINE init_symmetry(sym, atype, weight)
     944              :       TYPE(molsym_type), INTENT(inout)                   :: sym
     945              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: atype
     946              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: weight
     947              : 
     948              :       INTEGER                                            :: i, iatom, natoms
     949              : 
     950              :       ! Define the Cartesian axis vectors
     951          232 :       sym%x_axis = (/1.0_dp, 0.0_dp, 0.0_dp/)
     952          232 :       sym%y_axis = (/0.0_dp, 1.0_dp, 0.0_dp/)
     953          232 :       sym%z_axis = (/0.0_dp, 0.0_dp, 1.0_dp/)
     954              : 
     955              :       ! Initialize lists for symmetry elements
     956        93728 :       sym%sec(:, :, :) = 0.0_dp
     957       373288 :       sym%ses(:, :, :) = 0.0_dp
     958         4930 :       sym%sig(:, :) = 0.0_dp
     959              : 
     960          232 :       sym%center_of_mass(:) = 0.0_dp
     961              : 
     962              :       ! Initialize symmetry element counters
     963           58 :       sym%ncn = 1
     964           58 :       sym%nsn = 1
     965         1218 :       sym%nsec(:) = 0
     966         2378 :       sym%nses(:) = 0
     967           58 :       sym%nsig = 0
     968              : 
     969              :       ! Initialize logical variables
     970           58 :       sym%cubic = .FALSE.
     971           58 :       sym%dgroup = .FALSE.
     972           58 :       sym%igroup = .FALSE.
     973           58 :       sym%invers = .FALSE.
     974           58 :       sym%linear = .FALSE.
     975           58 :       sym%maxis = .FALSE.
     976           58 :       sym%ogroup = .FALSE.
     977           58 :       sym%sgroup = .FALSE.
     978           58 :       sym%sigmad = .FALSE.
     979           58 :       sym%sigmah = .FALSE.
     980           58 :       sym%sigmav = .FALSE.
     981           58 :       sym%tgroup = .FALSE.
     982              : 
     983              :       ! Initialize list of equivalent atoms (C1 symmetry)
     984           58 :       natoms = SIZE(sym%nequatom)
     985           58 :       sym%ngroup = natoms
     986          633 :       sym%nequatom(:) = 1
     987          633 :       DO i = 1, sym%ngroup
     988          575 :          sym%group_of(i) = i
     989          575 :          sym%llequatom(i) = i
     990          575 :          sym%symequ_list(i) = i
     991          633 :          sym%ulequatom(i) = i
     992              :       END DO
     993              : 
     994           58 :       sym%point_group_order = -1
     995              : 
     996          633 :       DO iatom = 1, natoms
     997          633 :          sym%aw(iatom) = weight(iatom)
     998              :       END DO
     999              : 
    1000              :       ! Generate atomic identification numbers (symmetry code) ***
    1001          633 :       sym%ain(:) = REAL(atype(:), KIND=dp) + sym%aw(:)
    1002              : 
    1003              :       ! Initialize the transformation matrix for input orientation -> standard orientation
    1004           58 :       CALL unit_matrix(sym%inptostd(:, :))
    1005              : 
    1006           58 :    END SUBROUTINE init_symmetry
    1007              : 
    1008              : ! **************************************************************************************************
    1009              : !> \brief  Return .TRUE., if the atom atoma is in the list symequ_list.
    1010              : !> \param atoma ...
    1011              : !> \param sym ...
    1012              : !> \return ...
    1013              : !> \par History
    1014              : !>        Creation (21.04.95, Matthias Krack)
    1015              : !> \author jgh
    1016              : ! **************************************************************************************************
    1017         1211 :    FUNCTION in_symequ_list(atoma, sym)
    1018              :       INTEGER, INTENT(IN)                                :: atoma
    1019              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1020              :       LOGICAL                                            :: in_symequ_list
    1021              : 
    1022              :       INTEGER                                            :: iatom
    1023              : 
    1024         1211 :       in_symequ_list = .FALSE.
    1025              : 
    1026         7969 :       DO iatom = 1, sym%ulequatom(sym%ngroup)
    1027         7969 :          IF (sym%symequ_list(iatom) == atoma) THEN
    1028         1211 :             in_symequ_list = .TRUE.
    1029              :             RETURN
    1030              :          END IF
    1031              :       END DO
    1032              : 
    1033              :    END FUNCTION in_symequ_list
    1034              : 
    1035              : ! **************************************************************************************************
    1036              : !> \brief Search for point groups of LOW SYMmetry (Ci,Cs).
    1037              : !> \param sym ...
    1038              : !> \param coord ...
    1039              : !> \par History
    1040              : !>        Creation (21.04.95, Matthias Krack)
    1041              : !> \author jgh
    1042              : ! **************************************************************************************************
    1043            3 :    SUBROUTINE lowsym(sym, coord)
    1044              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1045              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1046              : 
    1047              :       REAL(KIND=dp)                                      :: phi
    1048              :       REAL(KIND=dp), DIMENSION(3)                        :: d
    1049              : 
    1050            3 :       IF (sym%nsn == 2) THEN
    1051              :          ! Ci
    1052            1 :          sym%invers = .TRUE.
    1053            1 :          phi = angle(sym%ses(:, 1, 2), sym%z_axis(:))
    1054            1 :          d(:) = vector_product(sym%ses(:, 1, 2), sym%z_axis(:))
    1055            1 :          CALL rotate_molecule(phi, d(:), sym, coord)
    1056            2 :       ELSE IF (sym%nsig == 1) THEN
    1057              :          ! Cs
    1058            1 :          sym%sigmah = .TRUE.
    1059            1 :          phi = angle(sym%sig(:, 1), sym%z_axis(:))
    1060            1 :          d(:) = vector_product(sym%sig(:, 1), sym%z_axis(:))
    1061            1 :          CALL rotate_molecule(phi, d(:), sym, coord)
    1062              :       END IF
    1063              : 
    1064            3 :    END SUBROUTINE lowsym
    1065              : 
    1066              : ! **************************************************************************************************
    1067              : !> \brief Main program for symmetry analysis.
    1068              : !> \param sym ...
    1069              : !> \param eps_geo ...
    1070              : !> \param coord ...
    1071              : !> \param atype ...
    1072              : !> \param weight ...
    1073              : !> \par History
    1074              : !>        Creation (14.11.98, Matthias Krack)
    1075              : !> \author jgh
    1076              : ! **************************************************************************************************
    1077           58 :    SUBROUTINE molecular_symmetry(sym, eps_geo, coord, atype, weight)
    1078              :       TYPE(molsym_type), POINTER                         :: sym
    1079              :       REAL(KIND=dp), INTENT(IN)                          :: eps_geo
    1080              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1081              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: atype
    1082              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: weight
    1083              : 
    1084              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'molecular_symmetry'
    1085              : 
    1086              :       INTEGER                                            :: handle, natoms
    1087              : 
    1088           58 :       CALL timeset(routineN, handle)
    1089              : 
    1090              :       ! Perform memory allocation for the symmetry analysis
    1091           58 :       natoms = SIZE(coord, 2)
    1092           58 :       CALL create_molsym(sym, natoms)
    1093           58 :       sym%eps_geo = eps_geo
    1094              : 
    1095              :       ! Initialization of symmetry analysis
    1096           58 :       CALL init_symmetry(sym, atype, weight)
    1097              : 
    1098           58 :       IF (natoms == 1) THEN
    1099              :          ! Special case: only one atom
    1100            1 :          CALL atomsym(sym)
    1101              :       ELSE
    1102              :          ! Find molecular symmetry
    1103           57 :          CALL moleculesym(sym, coord)
    1104              :          ! Get point group and load point group table
    1105           57 :          CALL get_point_group_symbol(sym)
    1106              :       END IF
    1107              : 
    1108              :       ! Calculate group order
    1109           58 :       IF (.NOT. sym%linear) CALL get_point_group_order(sym)
    1110              : 
    1111              :       ! Generate a list of equivalent atoms
    1112           58 :       CALL build_symequ_list(sym, coord)
    1113              : 
    1114           58 :       CALL timestop(handle)
    1115              : 
    1116           58 :    END SUBROUTINE molecular_symmetry
    1117              : 
    1118              : ! **************************************************************************************************
    1119              : !> \brief Find the molecular symmetry.
    1120              : !> \param sym ...
    1121              : !> \param coord ...
    1122              : !> \par History
    1123              : !>        Creation (14.11.98, Matthias Krack)
    1124              : !> \author jgh
    1125              : ! **************************************************************************************************
    1126           57 :    SUBROUTINE moleculesym(sym, coord)
    1127              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1128              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1129              : 
    1130              :       INTEGER                                            :: icn, isec, isn
    1131              :       LOGICAL                                            :: failed
    1132              :       REAL(KIND=dp)                                      :: eps_tenval
    1133              : 
    1134              : ! Tolerance of degenerate eigenvalues for the molecular tensor of inertia
    1135              : 
    1136           57 :       eps_tenval = 0.01_dp*sym%eps_geo
    1137              : 
    1138              :       ! Calculate the molecular tensor of inertia
    1139           57 :       CALL tensor(sym, coord)
    1140              :       ! Use symmetry information from the eigenvalues of the molecular tensor of inertia
    1141           57 :       IF ((sym%tenval(3) - sym%tenval(1)) < eps_tenval) THEN
    1142              :          ! 0 < tenval(1) = tenval(2) = tenval(3)
    1143            7 :          sym%cubic = .TRUE.
    1144           50 :       ELSE IF ((sym%tenval(3) - sym%tenval(2)) < eps_tenval) THEN
    1145              :          ! 0 < tenval(1) < tenval(2) = tenval(3)
    1146              :          ! special case: 0 = tenval(1) < tenval(2) = tenval(3)
    1147           14 :          IF (sym%tenval(1) < eps_tenval) sym%linear = .TRUE.
    1148           14 :          CALL tracar(sym, coord, (/3, 1, 2/))
    1149              :       ELSE
    1150              :          ! 0 < tenval(1) = tenval(2) < tenval(3) or 0 < tenval(1) < tenval(2) < tenval(3)
    1151           36 :          CALL tracar(sym, coord, (/1, 2, 3/))
    1152              :       END IF
    1153              : 
    1154              :       ! Continue with the new coordinate axes
    1155              :       DO
    1156           57 :          failed = .FALSE.
    1157           57 :          IF (sym%cubic) THEN
    1158            7 :             CALL cubsym(sym, coord, failed)
    1159            7 :             IF (failed) THEN
    1160            0 :                sym%cubic = .FALSE.
    1161            0 :                CYCLE
    1162              :             END IF
    1163              : 
    1164           50 :          ELSE IF (sym%linear) THEN
    1165              : 
    1166              :             ! Linear molecule
    1167            2 :             IF (saxis(2, sym%z_axis(:), sym, coord)) THEN
    1168            1 :                sym%invers = .TRUE.
    1169            1 :                sym%dgroup = .TRUE.
    1170            1 :                sym%sigmah = .TRUE.
    1171              :             END IF
    1172              : 
    1173              :          ELSE
    1174              : 
    1175              :             ! Check the new coordinate axes for Cn axes
    1176          960 :             DO icn = 2, maxcn
    1177          912 :                IF (caxis(icn, sym%z_axis(:), sym, coord)) CALL addsec(icn, sym%z_axis(:), sym)
    1178          912 :                IF (caxis(icn, sym%x_axis(:), sym, coord)) CALL addsec(icn, sym%x_axis(:), sym)
    1179          960 :                IF (caxis(icn, sym%y_axis(:), sym, coord)) CALL addsec(icn, sym%y_axis(:), sym)
    1180              :             END DO
    1181              : 
    1182              :             ! Check the new coordinate axes for Sn axes
    1183         1920 :             DO isn = 2, maxsn
    1184         1872 :                IF (saxis(isn, sym%z_axis(:), sym, coord)) CALL addses(isn, sym%z_axis(:), sym)
    1185         1872 :                IF (saxis(isn, sym%x_axis(:), sym, coord)) CALL addses(isn, sym%x_axis(:), sym)
    1186         1920 :                IF (saxis(isn, sym%y_axis(:), sym, coord)) CALL addses(isn, sym%y_axis(:), sym)
    1187              :             END DO
    1188              : 
    1189              :             ! Check the new coordinate planes for mirror planes
    1190           48 :             IF (sigma(sym%z_axis(:), sym, coord)) CALL addsig(sym%z_axis(:), sym)
    1191           48 :             IF (sigma(sym%x_axis(:), sym, coord)) CALL addsig(sym%x_axis(:), sym)
    1192           48 :             IF (sigma(sym%y_axis(:), sym, coord)) CALL addsig(sym%y_axis(:), sym)
    1193              : 
    1194              :             ! There is a main axis (MAXIS = .TRUE.)
    1195           48 :             IF ((sym%ncn > 1) .OR. (sym%nsn > 3)) THEN
    1196           45 :                sym%maxis = .TRUE.
    1197           45 :                CALL axsym(coord, sym)
    1198              :             ELSE
    1199              :                ! Only low or no symmetry (Ci, Cs or C1)
    1200            3 :                CALL lowsym(sym, coord)
    1201              :             END IF
    1202              : 
    1203              :          END IF
    1204              : 
    1205           57 :          IF (.NOT. failed) EXIT
    1206              : 
    1207              :       END DO
    1208              : 
    1209              :       ! Find the remaining S axes
    1210           57 :       IF (.NOT. sym%linear) THEN
    1211          249 :          DO icn = 2, sym%ncn
    1212          545 :             DO isec = 1, sym%nsec(icn)
    1213          296 :                IF (saxis(icn, sym%sec(:, isec, icn), sym, coord)) &
    1214           89 :                   CALL addses(icn, sym%sec(:, isec, icn), sym)
    1215          296 :                IF (saxis(2*icn, sym%sec(:, isec, icn), sym, coord)) &
    1216          243 :                   CALL addses(2*icn, sym%sec(:, isec, icn), sym)
    1217              :             END DO
    1218              :          END DO
    1219              :       END IF
    1220              : 
    1221              :       ! Remove redundant S2 axes
    1222           57 :       IF (sym%nses(2) > 0) THEN
    1223           16 :          sym%nses(1) = sym%nses(1) - sym%nses(2)
    1224           16 :          sym%nses(2) = 0
    1225           16 :          CALL addses(2, sym%z_axis(:), sym)
    1226              :       END IF
    1227              : 
    1228           57 :    END SUBROUTINE moleculesym
    1229              : 
    1230              : ! **************************************************************************************************
    1231              : !> \brief The number of atoms which are placed on an axis defined by the vector a is returned.
    1232              : !> \param a ...
    1233              : !> \param coord ...
    1234              : !> \param sym ...
    1235              : !> \return ...
    1236              : !> \par History
    1237              : !>        Creation (16.11.98, Matthias Krack)
    1238              : !> \author jgh
    1239              : ! **************************************************************************************************
    1240           76 :    FUNCTION naxis(a, coord, sym)
    1241              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
    1242              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1243              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1244              :       INTEGER                                            :: naxis
    1245              : 
    1246              :       INTEGER                                            :: iatom, natoms
    1247              :       REAL(KIND=dp)                                      :: length_of_a, length_of_b, scapro
    1248              :       REAL(KIND=dp), DIMENSION(3)                        :: a_norm, b, b_norm
    1249              : 
    1250           76 :       naxis = 0
    1251           76 :       natoms = SIZE(coord, 2)
    1252              : 
    1253           76 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
    1254              : 
    1255              :       ! Check the length of vector a
    1256           76 :       IF (length_of_a > sym%eps_geo) THEN
    1257              : 
    1258          956 :          DO iatom = 1, natoms
    1259         3520 :             b(:) = coord(:, iatom)
    1260          880 :             length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
    1261              :             ! An atom in the origin counts for each axis
    1262          956 :             IF (length_of_b < sym%eps_geo) THEN
    1263            0 :                naxis = naxis + 1
    1264              :             ELSE
    1265         3520 :                a_norm = a(:)/length_of_a
    1266         3520 :                b_norm = b(:)/length_of_b
    1267          880 :                scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
    1268          880 :                IF (ABS(ABS(scapro) - 1.0_dp) < sym%eps_geo) naxis = naxis + 1
    1269              :             END IF
    1270              :          END DO
    1271              : 
    1272              :       END IF
    1273              : 
    1274           76 :    END FUNCTION naxis
    1275              : 
    1276              : ! **************************************************************************************************
    1277              : !> \brief Return .FALSE., if the symmetry elements se1 and se2 are a pair of
    1278              : !>         redundant symmetry elements.
    1279              : !> \param se1 ...
    1280              : !> \param se2 ...
    1281              : !> \return ...
    1282              : !> \par History
    1283              : !>        Creation (16.11.98, Matthias Krack)
    1284              : !> \author jgh
    1285              : ! **************************************************************************************************
    1286         1802 :    FUNCTION newse(se1, se2)
    1287              :       INTEGER, INTENT(IN)                                :: se1, se2
    1288              :       LOGICAL                                            :: newse
    1289              : 
    1290              :       INTEGER                                            :: ise
    1291              : 
    1292         1802 :       newse = .TRUE.
    1293              : 
    1294         3878 :       DO ise = 2, MIN(se1, se2)
    1295         3878 :          IF ((MODULO(se1, ise) == 0) .AND. (MODULO(se2, ise) == 0)) THEN
    1296         1802 :             newse = .FALSE.
    1297              :             RETURN
    1298              :          END IF
    1299              :       END DO
    1300              : 
    1301              :    END FUNCTION newse
    1302              : 
    1303              : ! **************************************************************************************************
    1304              : !> \brief The number of atoms which are placed in a plane defined by the
    1305              : !>         the normal vector a is returned.
    1306              : !> \param a ...
    1307              : !> \param sym ...
    1308              : !> \param coord ...
    1309              : !> \return ...
    1310              : !> \par History
    1311              : !>        Creation (16.11.98, Matthias Krack)
    1312              : !> \author jgh
    1313              : ! **************************************************************************************************
    1314          105 :    FUNCTION nsigma(a, sym, coord)
    1315              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
    1316              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1317              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1318              :       INTEGER                                            :: nsigma
    1319              : 
    1320              :       INTEGER                                            :: iatom, natoms
    1321              :       REAL(KIND=dp)                                      :: length_of_a, length_of_b, scapro
    1322              :       REAL(KIND=dp), DIMENSION(3)                        :: a_norm, b, b_norm
    1323              : 
    1324          105 :       nsigma = 0
    1325              : 
    1326          105 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
    1327              : 
    1328              :       ! Check the length of vector a
    1329          105 :       IF (length_of_a > sym%eps_geo) THEN
    1330          105 :          natoms = SIZE(coord, 2)
    1331          998 :          DO iatom = 1, natoms
    1332         3572 :             b(:) = coord(:, iatom)
    1333          893 :             length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
    1334              :             ! An atom in the origin counts for each mirror plane
    1335          998 :             IF (length_of_b < sym%eps_geo) THEN
    1336            0 :                nsigma = nsigma + 1
    1337              :             ELSE
    1338         3572 :                a_norm = a(:)/length_of_a
    1339         3572 :                b_norm = b(:)/length_of_b
    1340          893 :                scapro = a_norm(1)*b_norm(1) + a_norm(2)*b_norm(2) + a_norm(3)*b_norm(3)
    1341          893 :                IF (ABS(scapro) < sym%eps_geo) nsigma = nsigma + 1
    1342              :             END IF
    1343              :          END DO
    1344              :       END IF
    1345              : 
    1346          105 :    END FUNCTION nsigma
    1347              : 
    1348              : ! **************************************************************************************************
    1349              : !> \brief Style the output of the symmetry elements.
    1350              : !> \param se ...
    1351              : !> \param eps ...
    1352              : !> \par History
    1353              : !>        Creation (16.11.98, Matthias Krack)
    1354              : !> \author jgh
    1355              : ! **************************************************************************************************
    1356          522 :    SUBROUTINE outse(se, eps)
    1357              :       REAL(KIND=dp), DIMENSION(3), INTENT(INOUT)         :: se
    1358              :       REAL(KIND=dp), INTENT(IN)                          :: eps
    1359              : 
    1360              : ! Positive z component for all vectors
    1361              : 
    1362          522 :       IF (ABS(se(3)) < eps) THEN
    1363          255 :          IF (ABS(se(1)) < eps) THEN
    1364           47 :             se(2) = -se(2)
    1365              :          ELSE
    1366          448 :             IF (se(1) < 0.0_dp) se(1:2) = -se(1:2)
    1367              :          END IF
    1368              :       ELSE
    1369          471 :          IF (se(3) < 0.0_dp) se(:) = -se(:)
    1370              :       END IF
    1371              : 
    1372          522 :    END SUBROUTINE outse
    1373              : 
    1374              : ! **************************************************************************************************
    1375              : !> \brief Print the output of the symmetry analysis.
    1376              : !> \param sym ...
    1377              : !> \param coord ...
    1378              : !> \param atype ...
    1379              : !> \param element ...
    1380              : !> \param z ...
    1381              : !> \param weight ...
    1382              : !> \param iw ...
    1383              : !> \param plevel ...
    1384              : !> \par History
    1385              : !>        Creation (16.11.98, Matthias Krack)
    1386              : !> \author jgh
    1387              : ! **************************************************************************************************
    1388           58 :    SUBROUTINE print_symmetry(sym, coord, atype, element, z, weight, iw, plevel)
    1389              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1390              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(in)         :: coord
    1391              :       INTEGER, DIMENSION(:), INTENT(in)                  :: atype
    1392              :       CHARACTER(LEN=*), DIMENSION(:), INTENT(in)         :: element
    1393              :       INTEGER, DIMENSION(:), INTENT(in)                  :: z
    1394              :       REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: weight
    1395              :       INTEGER, INTENT(IN)                                :: iw, plevel
    1396              : 
    1397              :       REAL(KIND=dp), PARAMETER                           :: formatmaxval = 1.0E5_dp
    1398              : 
    1399              :       CHARACTER(LEN=3)                                   :: string
    1400              :       INTEGER                                            :: i, icn, iequatom, isec, ises, isig, isn, &
    1401              :                                                             j, nequatom, secount
    1402           58 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: equatomlist, idxlist
    1403              : 
    1404              :       ! Print point group symbol and point group order
    1405              :       WRITE (iw, "(/,(T2,A))") &
    1406           58 :          "MOLSYM| Molecular symmetry information", &
    1407          116 :          "MOLSYM|"
    1408              :       WRITE (iw, "(T2,A,T77,A4)") &
    1409           58 :          "MOLSYM| Point group", ADJUSTR(sym%point_group_symbol)
    1410           58 :       IF (sym%point_group_order > -1) THEN
    1411              :          WRITE (iw, "(T2,A,T77,I4)") &
    1412           55 :             "MOLSYM| Group order ", sym%point_group_order
    1413              :       END IF
    1414              : 
    1415           58 :       IF (MOD(plevel, 10) == 1) THEN
    1416              :          WRITE (iw, "(T2,A)") &
    1417           58 :             "MOLSYM|", &
    1418          116 :             "MOLSYM| Groups of atoms equivalent by symmetry"
    1419              :          WRITE (iw, "(T2,A,T31,A)") &
    1420           58 :             "MOLSYM| Group number", "Group members (atomic indices)"
    1421          164 :          DO i = 1, sym%ngroup
    1422          106 :             nequatom = sym%ulequatom(i) - sym%llequatom(i) + 1
    1423          106 :             CPASSERT(nequatom > 0)
    1424          424 :             ALLOCATE (equatomlist(nequatom), idxlist(nequatom))
    1425          681 :             equatomlist(1:nequatom) = sym%symequ_list(sym%llequatom(i):sym%ulequatom(i))
    1426          681 :             idxlist = 0
    1427          106 :             CALL sort(equatomlist, nequatom, idxlist)
    1428              :             WRITE (iw, "(T2,A,T10,I5,T16,I6,T31,10(1X,I4),/,(T2,'MOLSYM|',T31,10(1X,I4)))") &
    1429          106 :                "MOLSYM|", i, nequatom, (equatomlist(iequatom), iequatom=1, nequatom)
    1430          164 :             DEALLOCATE (equatomlist, idxlist)
    1431              :          END DO
    1432              :       END IF
    1433              : 
    1434           58 :       IF (MOD(plevel/100, 10) == 1) THEN
    1435              :          ! Print all symmetry elements
    1436              :          WRITE (iw, "(T2,A,/,T2,A,/,T2,A,T31,A,T45,A,T60,A,T75,A)") &
    1437           58 :             "MOLSYM|", &
    1438           58 :             "MOLSYM| Symmetry elements", &
    1439          116 :             "MOLSYM| Element number", "Type", "X", "Y", "Z"
    1440           58 :          secount = 1
    1441              :          WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A1,T36,3(1X,F14.8))") &
    1442           58 :             "MOLSYM|", secount, secount, "E", 0.0_dp, 0.0_dp, 0.0_dp
    1443              :          ! Mirror planes
    1444           58 :          string = "@  "
    1445          211 :          DO isig = 1, sym%nsig
    1446          153 :             secount = secount + 1
    1447          153 :             CALL outse(sym%sig(:, isig), sym%eps_geo)
    1448              :             WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
    1449          670 :                "MOLSYM|", secount, isig, string, (sym%sig(i, isig), i=1, 3)
    1450              :          END DO
    1451              :          ! C axes
    1452          252 :          DO icn = 2, sym%ncn
    1453          194 :             IF (icn < 10) THEN
    1454          194 :                WRITE (string, "(A1,I1)") "C", icn
    1455              :             ELSE
    1456            0 :                WRITE (string, "(A1,I2)") "C", icn
    1457              :             END IF
    1458          548 :             DO isec = 1, sym%nsec(icn)
    1459          296 :                secount = secount + 1
    1460          296 :                CALL outse(sym%sec(:, isec, icn), sym%eps_geo)
    1461              :                WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
    1462         1378 :                   "MOLSYM|", secount, isec, string, (sym%sec(i, isec, icn), i=1, 3)
    1463              :             END DO
    1464              :          END DO
    1465              :          ! S axes
    1466          215 :          DO isn = 2, sym%nsn
    1467          157 :             IF (isn == 2) THEN
    1468           29 :                WRITE (string, "(A3)") "i  "
    1469          128 :             ELSE IF (isn < 10) THEN
    1470          111 :                WRITE (string, "(A1,I1)") "S", isn
    1471              :             ELSE
    1472           17 :                WRITE (string, "(A1,I2)") "S", isn
    1473              :             END IF
    1474          288 :             DO ises = 1, sym%nses(isn)
    1475           73 :                secount = secount + 1
    1476           73 :                CALL outse(sym%ses(:, ises, isn), sym%eps_geo)
    1477              :                WRITE (iw, "(T2,A,T12,I5,T19,I5,T32,A3,T36,3(1X,F14.8))") &
    1478          449 :                   "MOLSYM|", secount, ises, string, (sym%ses(i, ises, icn), i=1, 3)
    1479              :             END DO
    1480              :          END DO
    1481              :       END IF
    1482              : 
    1483           58 :       IF (MOD(plevel/10, 10) == 1 .AND. SIZE(coord, 2) > 1) THEN
    1484              :          ! Print the moments of the molecular inertia tensor
    1485              :          WRITE (iw, "(T2,A,/,T2,A,/,T2,A,T43,A,T58,A,T73,A)") &
    1486           57 :             "MOLSYM|", &
    1487           57 :             "MOLSYM| Moments of the molecular inertia tensor [a.u.]", &
    1488          114 :             "MOLSYM|", "I(x)", "I(y)", "I(z)"
    1489           57 :          string = "xyz"
    1490          741 :          IF (MAXVAL(sym%tenmat) >= formatmaxval) THEN
    1491            0 :             DO i = 1, 3
    1492              :                WRITE (iw, "(T2,A,T32,A,T36,3(1X,ES14.4))") &
    1493            0 :                   "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
    1494              :             END DO
    1495              :          ELSE
    1496          228 :             DO i = 1, 3
    1497              :                WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
    1498          741 :                   "MOLSYM|", "I("//string(i:i)//")", (sym%tenmat(i, j), j=1, 3)
    1499              :             END DO
    1500              :          END IF
    1501              :          WRITE (iw, "(T2,A)") &
    1502           57 :             "MOLSYM|", &
    1503          114 :             "MOLSYM| Principal moments and axes of inertia [a.u.]"
    1504          741 :          IF (MAXVAL(sym%tenmat) >= formatmaxval) THEN
    1505              :             WRITE (iw, "(T2,A,T36,3(1X,ES14.4))") &
    1506            0 :                "MOLSYM|", (sym%tenval(i), i=1, 3)
    1507              :          ELSE
    1508              :             WRITE (iw, "(T2,A,T36,3(1X,F14.8))") &
    1509          228 :                "MOLSYM|", (sym%tenval(i), i=1, 3)
    1510              :          END IF
    1511          228 :          DO i = 1, 3
    1512              :             WRITE (iw, "(T2,A,T32,A,T36,3(1X,F14.8))") &
    1513          741 :                "MOLSYM|", string(i:i), (sym%tenvec(i, j), j=1, 3)
    1514              :          END DO
    1515              :       END IF
    1516              : 
    1517           58 :       IF (MOD(plevel, 10) == 1) THEN
    1518              :          ! Print the Cartesian coordinates of the standard orientation
    1519           58 :          CALL write_coordinates(coord, atype, element, z, weight, iw)
    1520              :       END IF
    1521              : 
    1522           58 :    END SUBROUTINE print_symmetry
    1523              : 
    1524              : ! **************************************************************************************************
    1525              : !> \brief Rotate the molecule about an axis defined by the vector a. The
    1526              : !>        rotation angle is phi (radians).
    1527              : !> \param phi ...
    1528              : !> \param a ...
    1529              : !> \param sym ...
    1530              : !> \param coord ...
    1531              : !> \par History
    1532              : !>        Creation (16.11.98, Matthias Krack)
    1533              : !> \author jgh
    1534              : ! **************************************************************************************************
    1535           87 :    SUBROUTINE rotate_molecule(phi, a, sym, coord)
    1536              :       REAL(KIND=dp), INTENT(IN)                          :: phi
    1537              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
    1538              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1539              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1540              : 
    1541              :       INTEGER                                            :: i
    1542              :       REAL(KIND=dp)                                      :: length_of_a
    1543              : 
    1544              :       ! Check the length of vector a
    1545           87 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
    1546           87 :       IF (length_of_a > sym%eps_geo) THEN
    1547              : 
    1548              :          ! Build up the rotation matrix
    1549           42 :          CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
    1550              : 
    1551              :          ! Rotate the molecule by phi around a
    1552        10689 :          coord(:, :) = MATMUL(sym%rotmat(:, :), coord(:, :))
    1553              : 
    1554              :          ! Rotate the normal vectors of the mirror planes which are already found
    1555         3339 :          sym%sig(:, 1:sym%nsig) = MATMUL(sym%rotmat(:, :), sym%sig(:, 1:sym%nsig))
    1556              : 
    1557              :          ! Rotate the Cn axes which are already found
    1558          186 :          DO i = 2, sym%ncn
    1559         6843 :             sym%sec(:, 1:sym%nsec(i), i) = MATMUL(sym%rotmat(:, :), sym%sec(:, 1:sym%nsec(i), i))
    1560              :          END DO
    1561              : 
    1562              :          ! Rotate the Sn axes which are already found
    1563          134 :          DO i = 2, sym%nsn
    1564         2150 :             sym%ses(:, 1:sym%nses(i), i) = MATMUL(sym%rotmat(:, :), sym%ses(:, 1:sym%nses(i), i))
    1565              :          END DO
    1566              : 
    1567              :          ! Store current rotation
    1568         2688 :          sym%inptostd(:, :) = MATMUL(sym%rotmat(:, :), sym%inptostd(:, :))
    1569              : 
    1570              :       END IF
    1571              : 
    1572           87 :    END SUBROUTINE rotate_molecule
    1573              : 
    1574              : ! **************************************************************************************************
    1575              : !> \brief Rotate the molecule about a n-fold axis defined by the vector a
    1576              : !>        using a rotation angle of 2*pi/n. Check, if the axis is a Cn axis.
    1577              : !> \param n ...
    1578              : !> \param a ...
    1579              : !> \param sym ...
    1580              : !> \param coord ...
    1581              : !> \return ...
    1582              : !> \par History
    1583              : !>        Creation (16.11.98, Matthias Krack)
    1584              : !> \author jgh
    1585              : ! **************************************************************************************************
    1586         6344 :    FUNCTION saxis(n, a, sym, coord)
    1587              :       INTEGER, INTENT(IN)                                :: n
    1588              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
    1589              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1590              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1591              :       LOGICAL                                            :: saxis
    1592              : 
    1593              :       INTEGER                                            :: iatom, natoms
    1594              :       REAL(KIND=dp)                                      :: length_of_a, phi
    1595              :       REAL(KIND=dp), DIMENSION(3)                        :: b
    1596              : 
    1597         6344 :       saxis = .FALSE.
    1598              : 
    1599         6344 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
    1600              : 
    1601         6344 :       natoms = SIZE(coord, 2)
    1602              : 
    1603              :       ! Check the length of the axis vector a
    1604         6344 :       IF (length_of_a > sym%eps_geo) THEN
    1605              :          ! Calculate the rotation angle phi and build up the rotation matrix rotmat
    1606         6344 :          phi = 2.0_dp*pi/REAL(n, KIND=dp)
    1607         6344 :          CALL build_rotmat(phi, a(:), sym%rotmat(:, :))
    1608              :          ! Check all atoms
    1609         9816 :          DO iatom = 1, natoms
    1610              :             ! Rotate the actual atom by 2*pi/n about a
    1611       124111 :             b(:) = MATMUL(sym%rotmat(:, :), coord(:, iatom))
    1612              :             ! Reflect the coordinates of the rotated atom
    1613        38188 :             b(:) = reflect_vector(b(:), a(:))
    1614              :             ! Check if the coordinates of the rotated atom are in the coordinate set of the molecule
    1615         9816 :             IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
    1616              :          END DO
    1617              :          ! The axis defined by a is a Sn axis, thus return with saxis = .TRUE.
    1618              :          saxis = .TRUE.
    1619              :       END IF
    1620              : 
    1621              :    END FUNCTION saxis
    1622              : 
    1623              : ! **************************************************************************************************
    1624              : !> \brief Reflect all atoms of the molecule through a mirror plane defined
    1625              : !>        by the normal vector a.
    1626              : !> \param a ...
    1627              : !> \param sym ...
    1628              : !> \param coord ...
    1629              : !> \return ...
    1630              : !> \par History
    1631              : !>        Creation (16.11.98, Matthias Krack)
    1632              : !> \author jgh
    1633              : ! **************************************************************************************************
    1634         2808 :    FUNCTION sigma(a, sym, coord)
    1635              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a
    1636              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1637              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1638              :       LOGICAL                                            :: sigma
    1639              : 
    1640              :       INTEGER                                            :: iatom, natoms
    1641              :       REAL(KIND=dp)                                      :: length_of_a
    1642              :       REAL(KIND=dp), DIMENSION(3)                        :: b
    1643              : 
    1644         2808 :       sigma = .FALSE.
    1645              : 
    1646         2808 :       length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
    1647              : 
    1648              :       ! Check the length of vector a
    1649         2808 :       IF (length_of_a > sym%eps_geo) THEN
    1650         2674 :          natoms = SIZE(coord, 2)
    1651        10068 :          DO iatom = 1, natoms
    1652              :             ! Reflect the actual atom
    1653         9471 :             b(:) = reflect_vector(coord(:, iatom), a(:))
    1654              :             ! Check if the coordinates of the reflected atom are in the coordinate set of the molecule
    1655        10068 :             IF (equatom(iatom, b(:), sym, coord) == 0) RETURN
    1656              :          END DO
    1657              :          ! The plane defined by a is a mirror plane, thus return with sigma = .TRUE.
    1658              :          sigma = .TRUE.
    1659              :       END IF
    1660              : 
    1661              :    END FUNCTION sigma
    1662              : 
    1663              : ! **************************************************************************************************
    1664              : !> \brief Calculate the molecular tensor of inertia and the vector to the
    1665              : !>        center of mass of the molecule.
    1666              : !> \param sym ...
    1667              : !> \param coord ...
    1668              : !> \par History
    1669              : !>        Creation (16.11.98, Matthias Krack)
    1670              : !> \author jgh
    1671              : ! **************************************************************************************************
    1672           57 :    SUBROUTINE tensor(sym, coord)
    1673              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1674              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1675              : 
    1676              :       INTEGER                                            :: natoms
    1677              :       REAL(KIND=dp)                                      :: tt
    1678              : 
    1679              :       ! Find the vector center_of_mass to molecular center of mass
    1680           57 :       natoms = SIZE(coord, 2)
    1681         3269 :       sym%center_of_mass(:) = MATMUL(coord(:, 1:natoms), sym%aw(1:natoms))/SUM(sym%aw(1:natoms))
    1682              : 
    1683              :       ! Translate the center of mass of the molecule to the origin
    1684         2353 :       coord(:, 1:natoms) = coord(:, 1:natoms) - SPREAD(sym%center_of_mass(:), 2, natoms)
    1685              : 
    1686              :       ! Build up the molecular tensor of inertia
    1687              : 
    1688          631 :       sym%tenmat(1, 1) = DOT_PRODUCT(sym%aw(1:natoms), (coord(2, 1:natoms)**2 + coord(3, 1:natoms)**2))
    1689          631 :       sym%tenmat(2, 2) = DOT_PRODUCT(sym%aw(1:natoms), (coord(3, 1:natoms)**2 + coord(1, 1:natoms)**2))
    1690          631 :       sym%tenmat(3, 3) = DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)**2 + coord(2, 1:natoms)**2))
    1691              : 
    1692          631 :       sym%tenmat(1, 2) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(2, 1:natoms)))
    1693          631 :       sym%tenmat(1, 3) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(1, 1:natoms)*coord(3, 1:natoms)))
    1694          631 :       sym%tenmat(2, 3) = -DOT_PRODUCT(sym%aw(1:natoms), (coord(2, 1:natoms)*coord(3, 1:natoms)))
    1695              : 
    1696              :       ! Symmetrize tensor matrix
    1697           57 :       sym%tenmat(2, 1) = sym%tenmat(1, 2)
    1698           57 :       sym%tenmat(3, 1) = sym%tenmat(1, 3)
    1699           57 :       sym%tenmat(3, 2) = sym%tenmat(2, 3)
    1700              : 
    1701              :       ! Diagonalize the molecular tensor of inertia
    1702           57 :       CALL jacobi(sym%tenmat(:, :), sym%tenval(:), sym%tenvec(:, :))
    1703              : 
    1704              :       ! Secure that the principal axes are right-handed
    1705          228 :       sym%tenvec(:, 3) = vector_product(sym%tenvec(:, 1), sym%tenvec(:, 2))
    1706              : 
    1707           57 :       tt = SQRT(sym%tenval(1)**2 + sym%tenval(2)**2 + sym%tenval(3)**2)
    1708           57 :       CPASSERT(tt /= 0)
    1709              : 
    1710           57 :    END SUBROUTINE tensor
    1711              : 
    1712              : ! **************************************************************************************************
    1713              : !> \brief Transformation the Cartesian coodinates with the matrix tenvec.
    1714              : !>        The coordinate axes can be switched according to the index
    1715              : !>        vector idx. If idx(1) is equal to 3 instead to 1, then the first
    1716              : !>        eigenvector (smallest eigenvalue) of the molecular tensor of
    1717              : !>        inertia is connected to the z axis instead of the x axis. In
    1718              : !>        addition the local atomic coordinate systems are initialized,
    1719              : !>        if the symmetry is turned off.
    1720              : !> \param sym ...
    1721              : !> \param coord ...
    1722              : !> \param idx ...
    1723              : !> \par History
    1724              : !>        Creation (16.11.98, Matthias Krack)
    1725              : !> \author jgh
    1726              : ! **************************************************************************************************
    1727           50 :    SUBROUTINE tracar(sym, coord, idx)
    1728              :       TYPE(molsym_type), INTENT(inout)                   :: sym
    1729              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(inout)      :: coord
    1730              :       INTEGER, DIMENSION(3), INTENT(IN)                  :: idx
    1731              : 
    1732              :       INTEGER                                            :: iatom, natom
    1733              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: tenvec
    1734              : 
    1735              :       ! Transformation of the atomic coordinates  ***
    1736           50 :       natom = SIZE(coord, 2)
    1737          650 :       tenvec = TRANSPOSE(sym%tenvec)
    1738          482 :       DO iatom = 1, natom
    1739         3074 :          coord(idx, iatom) = MATMUL(tenvec, coord(:, iatom))
    1740              :       END DO
    1741              : 
    1742              :       ! Modify the transformation matrix for input orientation -> standard orientation
    1743          650 :       sym%inptostd(idx, :) = tenvec
    1744              : 
    1745           50 :    END SUBROUTINE tracar
    1746              : 
    1747              : ! **************************************************************************************************
    1748              : !> \brief Distance between two points
    1749              : !> \param a ...
    1750              : !> \param b ...
    1751              : !> \return ...
    1752              : !> \author jgh
    1753              : ! **************************************************************************************************
    1754         2864 :    FUNCTION dist(a, b) RESULT(d)
    1755              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: a, b
    1756              :       REAL(KIND=dp)                                      :: d
    1757              : 
    1758        11456 :       d = SQRT(SUM((a - b)**2))
    1759              : 
    1760         2864 :    END FUNCTION
    1761              : ! **************************************************************************************************
    1762              : !> \brief   Write atomic coordinates to output unit iw.
    1763              : !> \param coord ...
    1764              : !> \param atype ...
    1765              : !> \param element ...
    1766              : !> \param z ...
    1767              : !> \param weight ...
    1768              : !> \param iw ...
    1769              : !> \date    08.05.2008
    1770              : !> \author  JGH
    1771              : ! **************************************************************************************************
    1772           58 :    SUBROUTINE write_coordinates(coord, atype, element, z, weight, iw)
    1773              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(in)         :: coord
    1774              :       INTEGER, DIMENSION(:), INTENT(in)                  :: atype
    1775              :       CHARACTER(LEN=*), DIMENSION(:), INTENT(in)         :: element
    1776              :       INTEGER, DIMENSION(:), INTENT(in)                  :: z
    1777              :       REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: weight
    1778              :       INTEGER, INTENT(IN)                                :: iw
    1779              : 
    1780              :       INTEGER                                            :: iatom, natom
    1781              : 
    1782           58 :       IF (iw > 0) THEN
    1783           58 :          natom = SIZE(coord, 2)
    1784              :          WRITE (UNIT=iw, FMT="(T2,A)") &
    1785           58 :             "MOLSYM|", &
    1786          116 :             "MOLSYM| Atomic coordinates of the standard orientation in BOHR"
    1787              :          WRITE (UNIT=iw, FMT="(T2,A,T37,A,T51,A,T65,A,T77,A)") &
    1788           58 :             "MOLSYM|   Atom Kind Element", "X", "Y", "Z", "Mass"
    1789          633 :          DO iatom = 1, natom
    1790              :             WRITE (UNIT=iw, FMT="(T2,A,I7,1X,I4,1X,A2,1X,I4,3(1X,F13.8),1X,F9.4)") &
    1791          575 :                "MOLSYM|", iatom, atype(iatom), ADJUSTL(element(iatom)), z(iatom), &
    1792         1208 :                coord(1:3, iatom), weight(iatom)
    1793              :          END DO
    1794              :          WRITE (UNIT=iw, FMT="(T2,A)") &
    1795           58 :             "MOLSYM|", &
    1796          116 :             "MOLSYM| Atomic coordinates of the standard orientation in ANGSTROM"
    1797              :          WRITE (UNIT=iw, FMT="(T2,A,T37,A,T51,A,T65,A,T77,A)") &
    1798           58 :             "MOLSYM|   Atom Kind Element", "X", "Y", "Z", "Mass"
    1799          633 :          DO iatom = 1, natom
    1800              :             WRITE (UNIT=iw, FMT="(T2,A,I7,1X,I4,1X,A2,1X,I4,3(1X,F13.8),1X,F9.4)") &
    1801          575 :                "MOLSYM|", iatom, atype(iatom), ADJUSTL(element(iatom)), z(iatom), &
    1802         2933 :                coord(1:3, iatom)*angstrom, weight(iatom)
    1803              :          END DO
    1804              :       END IF
    1805              : 
    1806           58 :    END SUBROUTINE write_coordinates
    1807              : 
    1808            0 : END MODULE molsym
        

Generated by: LCOV version 2.0-1