C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%   Subroutine "symmetry" is used to get the symmetry      %%%%
C%%%%   (i.e. the space group number) from the structrue file  %%%%
C%%%%   using the spglib, which is a library for finding and   %%%%
C%%%%   handling crystal symmetries                            %%%%
C%%%%                           REFERENCE                      %%%%
C%%%%   1\ https://atztogo.github.io/spglib/                   %%%% 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        module defs_basis
          implicit none
          integer, parameter :: dp=kind(1.0d0)
        end module defs_basis

        subroutine symmetry(spgnum, inpos)

        use defs_basis
        use spglib_f08
        use constant
        implicit none        

        integer spgnum

        double precision alat,pripos,privect
        dimension privect(3,3),pripos(1000,3)
        character*10 title,fix
        character*1 pricoord,seldyn
        character*5 sname
        character*10 inpos   
        dimension sname(10),fix(1000,3)
        integer i,j,k,ntype,natomi,nn
        integer Ncoord,Ndyn,atom_types,anum
        dimension natomi(10),atom_types(1000),anum(11) 

        character*11 symbol
        integer space_group
        real*8 lattice, positions
        dimension lattice(3,3), positions(3,1000)


C       >>> Tolerance of distance between atomic positions and 
C           between lengths of lattice vectors to be tolerated 
C           in the symmetry finding.      

        symprec=0.01        

C       >>> End

        call readpos(inpos,
     &               title, 
     &               alat,
     &               privect, 
     &               sname,
     &               natomi,                
     &               ntype,
     &               nn,  
     &               Ndyn,
     &               Ncoord,
     &               pripos,
     &               fix,
     &               atom_types)


        if (Ncoord .EQ. 1 .OR. Ncoord .EQ. 0) THEN
                call directpos(alat,
     &                        privect,
     &                        nn,
     &                        Ncoord,
     &                        pripos)
                Ncoord=0
        else
                open(unit=19,position='Append',FILE='RELAS')
                write(19,*) 'ERROR!! NOT DIRECT OR CARTESIAN!!!'
                close(19)
                stop
        endif

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%    To get the space group of the structure "INPOS"      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        do i=1,3
           do j=1,3
              lattice(i,j)=privect(i,j)
           enddo
        enddo

        do i=1,nn
           do j=1,3
              positions(j,i)=pripos(i,j)
           enddo
        enddo

        symbol='             '
        space_group = spg_get_international( symbol, lattice, 
     &                positions, atom_types, nn, symprec )

c        open(unit=19,position='Append',FILE='RELAS') 
        
C        write(19,'(a28,a,a1,i03.3,a1)') 'Space group symbol(number): '
C     &                           , trim(symbol), '(', space_group, ')'

C        close(19)

        if ( space_group .EQ. 0) then
                open(unit=19,position='Append',FILE='RELAS') 
                write(19,*) 'ERROR!! Cannot get space group!!!'
                close(19)
                stop

        else
                spgnum=space_group
        endif
        
        end 
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                    End of symmetry.f                    %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
