program detsym
      implicit none
      real*8,allocatable::r(:),line(:),u(:,:)
      integer*4,allocatable::q(:),table(:,:)
      integer*4 ns0,i,nat,n,j,ix,iy,iz,mx,n6,ia,iap,ii,ee,ib,ie,ifa,
     1ig,im,ixa,nb,nt,coord,ixp,iip,jj
      parameter (ns0=10,mx=12,n6=6)
      real*8 o(ns0,3,3),x,y,z,xs,ys,zs,xp,yp,zp,tol,d,ln,sn
      character*3 sst(ns0),SSP(ns0,mx),so(ns0)
      character*1 xyz(3)
      data xyz/'x','y','z'/
      integer*4 CHT(ns0,mx,mx),KG(ns0),se(ns0,mx)
      logical debug

      write(6,598)
 598  format(/,' Symmetry determination',/)

c     distance limit = 10-2A, tol=(10-2)^2
      tol=1.0d-4
      debug=.false.

      open(9,file='FILE.X')
      read(9,*)
      read(9,*)nat
      n=3*nat
      allocate(r(n),q(nat))
      do 1 i=1,nat
  1   read(9,*)q(i),(r(j+3*(i-1)),j=1,3)
      close(9)
      write(6,599)nat
 599  format(i5,' coord atoms in FILE.X',/)

      call setel(ns0,o,so)
      call SETSYM(CHT,sst,KG,SSP,se)

c     loop over symmetry point groups:
c     find one with maximum symmetry elements:
      ifa=0
      im=0
      do 21 ig=1,n6

      if(debug)write(6,6001)ig,sst(ig)
6001  format(i2,' group ',A3) 

c     try all axis combinations:
      do 21 ix=1,3
      iy=ix+1
      if(iy.eq.4)iy=1
      iz=iy+1
      if(iz.eq.4)iz=1

      if(debug)write(6,6002)ix,iy,iz
6002  format('xyz:',3i2)

c     loop over symmetry elements:
      do 2 ie=1,KG(ig)
      ee=se(ig,ie)
      if(debug)write(6,6003)ie,ee
6003  format(i2,' element,',i2,' operation')
c     loop over atoms:
      do 31 ia=1,nat
      x=r(ix+3*(ia-1))
      y=r(iy+3*(ia-1))
      z=r(iz+3*(ia-1))
c     symmetry-related coordinate
      xs=o(ee,1,1)*x+o(ee,1,2)*y+o(ee,1,3)*z
      ys=o(ee,2,1)*x+o(ee,2,2)*y+o(ee,2,3)*z
      zs=o(ee,3,1)*x+o(ee,3,2)*y+o(ee,3,3)*z
      if(debug)write(6,6004)ia,x,y,z,xs,ys,zs
6004  format(i2,'atom',6f10.3)

c     is in this position of an atom?:
      do 3 iap=1,nat
      xp=r(ix+3*(iap-1))
      yp=r(iy+3*(iap-1))
      zp=r(iz+3*(iap-1))
      d=(xp-xs)**2+(yp-ys)**2+(zp-zs)**2
c     it is, go to another atom: 
3     if(d.le.tol.and.q(ia).eq.q(iap))goto 31
      if(debug)write(6,*)'not found'
c     it is not, try another symmetry group:
      goto 21
31    continue
c     ia
2     continue
c     ie

c     now for all symmetry elements atoms could be assigned
      write(6,601)sst(ig),KG(ig)
601   format(' ',a3,' group found, Nel = ',i2)
      if(KG(ig).gt.im)then
       im=KG(ig)
       ifa=ig
       ixa=ix
      endif
21    continue
c     ig,ix


      ix=ixa
      iy=ix+1
      if(iy.eq.4)iy=1
      iz=iy+1
      if(iz.eq.4)iz=1
      write(6,602)sst(ifa)
602   format(/,' ',a3,' is the best')
c     where atoms go under symmetry operations:
      allocate(table(n,im))
      if(debug)write(6,6002)ix,iy,iz
      table=0
      do 4 ie=1,KG(ifa)
      ee=se(ifa,ie)
      if(debug)write(6,6003)ie,ee
c     loop over atoms:
      do 4 ia=1,nat
      ii=3*(ia-1)
      x=r(ix+ii)
      y=r(iy+ii)
      z=r(iz+ii)
c     symmetry-related coordinate
      xs=o(ee,1,1)*x+o(ee,1,2)*y+o(ee,1,3)*z
      ys=o(ee,2,1)*x+o(ee,2,2)*y+o(ee,2,3)*z
      zs=o(ee,3,1)*x+o(ee,3,2)*y+o(ee,3,3)*z
      if(debug)write(6,6004)ia,x,y,z,xs,ys,zs
      do 4 iap=1,nat
      xp=r(ix+3*(iap-1))
      yp=r(iy+3*(iap-1))
      zp=r(iz+3*(iap-1))
      d=(xs-xp)**2+(ys-yp)**2+(zs-zp)**2
      if(d.le.tol.and.q(ia).eq.q(iap))then
c      look at symmetry-related coordinate changes
       if(o(ee,1,1).gt.0.0d0)then
        table(ix+ii,ie)= ix+3*(iap-1)
       else
        table(ix+ii,ie)=-ix-3*(iap-1)
       endif
       if(o(ee,2,2).gt.0.0d0)then
        table(iy+ii,ie)= iy+3*(iap-1)
       else
        table(iy+ii,ie)=-iy-3*(iap-1)
       endif
       if(o(ee,3,3).gt.0.0d0)then
        table(iz+ii,ie)= iz+3*(iap-1)
       else
        table(iz+ii,ie)=-iz-3*(iap-1)
       endif
      endif
4     continue

      write(6,603)(so(se(ifa,i)),i=1,KG(ifa))
603   format(/,' coord atom xyz  ',12(a3,' '))
      do 5 ia=1,nat
      do 5 ix=1,3
      ii=ix+3*(ia-1)
5     write(6,604)ii,ia,xyz(ix),(table(ii,i),i=1,KG(ifa))
604   format(2i7,3x,a1,2x,12i6)

c     transformation matrix:
      allocate(u(n,n))
      u=0.0d0
c     symmetry blocks:
      allocate(line(n))
      nt=0
      do 6 ib=1,KG(ifa)
      write(6,605)SSP(ifa,ib)
605   format(' ',a3)
      nb=0
      do 71 ia=1,nat
      do 71 ix=1,3
      ii=ix+3*(ia-1)
      line=0.0d0
      do 7 i=1,KG(ifa)
      coord=iabs(table(ii,i))
      if(table(ii,i).lt.0)then
       sn=-1.0d0
      else
       sn= 1.0d0
      endif
7     line(coord)=line(coord)+sn*dble(CHT(ifa,ib,i))
      ln=0.0d0
      do 8 ii=1,n
 8    ln=ln+line(ii)**2
      ln=dsqrt(ln)
      if(ln.gt.tol)then
       do 9 ii=1,n
 9     line(ii)=line(ii)/ln


       nb=nb+1
       nt=nt+1
       write(6,606)nb
606    format(i3,$)
       do 10 iap=1,nat
       do 10 ixp=1,3
       ii=ixp+3*(iap-1)
 10    if(dabs(line(ii)).gt.tol)write(6,607)line(ii),iap,xyz(ixp)
607    format(f10.3,i6,a1,$)

c      check if already defined
       do 11 iip=1,nt-1
       sn=0.0d0
       do 111 jj=1,n
111    sn=sn+line(jj)*u(iip,jj)
       if(dabs(sn).gt.tol)then
        write(6,610)iip
610     format(' same as number ',i3,', deleted')
        nb=nb-1
        nt=nt-1
        goto 71
       endif
11     continue

       if(nt.le.n)then
        do 12 jj=1,n
12      u(nt,jj)=line(jj)
        write(6,*)' added'
       else
        write(6,*)'Error: too many symmetrized coordinates'
        stop
       endif
      endif
71    continue
c     ia,ix
6     write(6,608)nb
c     ib

608   format(i4,' elements')
      write(6,*)
      write(6,609)nt
609   format(i4,' in total')
      end

      SUBROUTINE SETSYM(CHT,sst,KG,SSP,se)
      integer*4 ns0,mx
      parameter (ns0=10,mx=12)
      integer*4 CHT(ns0,mx,mx),KG(ns0),se(ns0,mx)
      character*3 sst(ns0),SSP(ns0,mx)
      character*1 PRIME
      PRIME=CHAR(39)                               
      CHT=1
      sst='xxx'
      se=0
c     sst .. point group
c     KG .. number of symmetry elements
c     SSP .. symmetry

      sst(1)='C1 '
      KG(1)=1                 
      SSP(1,1)='A  '
      se(1,1)=1

      sst(2)='C2V'
      KG(2)=4                 
      CHT(2,2,3)=-1                                  
      CHT(2,2,4)=-1                                  
      CHT(2,3,2)=-1                                  
      CHT(2,3,4)=-1                                  
      CHT(2,4,2)=-1                                  
      CHT(2,4,3)=-1                                  
C       1 1 1 1 A1                                 
C       1 1-1-1 A2                                 
C       1-1 1-1 B1                                 
C       1-1-1 1 B2                                 
      SSP(2,1)='A1 '  
      SSP(2,2)='A2 ' 
      SSP(2,3)='B1 ' 
      SSP(2,4)='B2 ' 
c     E C2Z sigmaxz sigmazy
      se(2,1)=1
      se(2,2)=2
      se(2,3)=3
      se(2,4)=4

      sst(3)='C2H'
      KG(3)=4                 
      CHT(3,2,2)=-1                                  
      CHT(3,2,4)=-1                                  
      CHT(3,3,3)=-1                                  
      CHT(3,3,4)=-1                                  
      CHT(3,4,2)=-1                                  
      CHT(3,4,3)=-1                                  
C       1 1 1 1 Ag                                 
C       1-1 1-1 Bg                                 
C       1 1-1-1 Au                                 
C       1-1-1 1 Bu                                 
      SSP(3,1)='Ag '                               
      SSP(3,2)='Bg '                               
      SSP(3,3)='Au '                               
      SSP(3,4)='Bu '                               
c     E C2Z i sigmah
      se(3,1)=1
      se(3,2)=2
      se(3,3)=5
      se(3,4)=6

      sst(4)='C2 '
      KG(4)=2                                        
      CHT(4,2,2)=-1                                 
C       1 1 A                                      
C       1-1 B                                      
      SSP(4,1)='A  '                                
      SSP(4,2)='B  '                                
c     E C2
      se(4,1)=1
      se(4,2)=2

      sst(5)='CS '
      KG(5)=2                                        
      CHT(5,2,2)=-1                                 
C       1 1 A'                                     
C       1-1 A''                                    
      SSP(5,1)='AP '                                
      SSP(5,1)(2:2)=PRIME                           
      SSP(5,2)='APP'                                
      SSP(5,2)(2:2)=PRIME                           
      SSP(5,2)(3:3)=PRIME                           
c     E sigmah
      se(5,1)=1
      se(5,2)=6
     
      sst(6)='D2H'
      KG(6)=8
      CHT(6,2,3)=-1                                 
      CHT(6,2,4)=-1                                 
      CHT(6,2,7)=-1                                 
      CHT(6,2,8)=-1                                 
      CHT(6,3,2)=-1                                 
      CHT(6,3,4)=-1                                 
      CHT(6,3,6)=-1                                 
      CHT(6,3,8)=-1                                 
      CHT(6,4,2)=-1                                 
      CHT(6,4,3)=-1                                 
      CHT(6,4,6)=-1                                 
      CHT(6,4,7)=-1                                 
      CHT(6,5,5)=-1                                 
      CHT(6,5,6)=-1                                 
      CHT(6,5,7)=-1                                 
      CHT(6,5,8)=-1                                 
      CHT(6,6,3)=-1                                 
      CHT(6,6,4)=-1                                 
      CHT(6,6,5)=-1                                 
      CHT(6,6,6)=-1                                 
      CHT(6,7,2)=-1                                 
      CHT(6,7,4)=-1                                 
      CHT(6,7,5)=-1                                 
      CHT(6,7,7)=-1                                 
      CHT(6,8,2)=-1                                 
      CHT(6,8,3)=-1                                 
      CHT(6,8,5)=-1                                 
      CHT(6,8,8)=-1                                 
C      1 1 1 1 1 1 1 1 Ag   1                      
C      1 1-1-1 1 1-1-1 B1g  2                      
C      1-1 1-1 1-1 1-1 B2g  3                      
C      1-1-1 1 1-1-1 1 B3g  4                      
C      1 1 1 1-1-1-1-1 Au   5                      
C      1 1-1-1-1-1 1 1 B1u  6                      
C      1-1 1-1-1 1-1 1 B2u  7                      
C      1-1-1 1-1 1 1-1 B3u  8                      
      SSP(6,1)='Ag '                                
      SSP(6,2)='B1g'                                
      SSP(6,3)='B2g'                                
      SSP(6,4)='B3g'                                
      SSP(6,5)='Au '                                
      SSP(6,6)='B1u'                                
      SSP(6,7)='B2u'                                
      SSP(6,8)='B3u'                                
c     E c2z c2y c2x i xy xz yz
      se(6,1)=1
      se(6,2)=2
      se(6,3)=7
      se(6,4)=8
      se(6,5)=5
      se(6,6)=6
      se(6,7)=3
      se(6,8)=4
      RETURN
      END  
      
      subroutine setel(ns0,o,so)
      implicit none
      integer*4 ns0,j,i
      real*8 o(ns0,3,3)
      character*3 so(ns0)
      o=0.0d0
c     1 identity, E:
      do 1 i=1,ns0
      do 1 j=1,3
1     o(i,j,j) =  1.0d0
      so(1)='E  '
c     2 C2z:
      o(2,1,1) = -1.0d0
      o(2,2,2) = -1.0d0
      so(2)='C2z'
c     3 sigmazx:
      o(3,2,2) = -1.0d0
      so(3)='szx'
c     4 sigmazy:
      o(4,1,1) = -1.0d0
      so(4)='szy'
c     5 inversion
      o(5,1,1) = -1.0d0
      o(5,2,2) = -1.0d0
      o(5,3,3) = -1.0d0
      so(5)='i  '
c     6 sigmaxy:
      o(6,3,3) = -1.0d0
      so(6)='sxy'
c     7 C2y:
      o(7,1,1) = -1.0d0
      o(7,3,3) = -1.0d0
      so(7)='C2y'
c     8 C2x:
      o(8,2,2) = -1.0d0
      o(8,3,3) = -1.0d0
      so(8)='C2x'
      return
      end