program read_incommensurate_subgroups
* read the data file containing the isotropy subgroups
      implicit none
      integer i,j,k,m,n,ng
      real sign
* space group number (1-230)
      integer ispacegroup(7799)
* space group symbol: P1, P-1, P2, etc.
      character spacegroupstring(7799)*10
* IR number (1-5508)
      integer irnumber(7799)
* IR symbol: GM1, DT1, etc.
      character irstring(7799)*8
* k vector symbol: GM, DT, etc.
      character kstring(7799)*2
* k vector string:  (a,b,g), (1/2,0,g), etc.
      character kvectorstring(7799)*12
* k vector components:  
* kvector(i,1,n) is the ith rational component
* kvector(i,2,n) is the ith alpha component
* kvector(i,3,n) is the ith beta component
* kvector(i,4,n) is the ith gamma component
      real kvector(3,4,7799)
* order parameter direction symbol
      character opstring(7799)*8
* order parameter direction string
      character opvectorstring(7799)*19
* dimension of IR
      integer irdim(7799)
* subduction frequency: number of degrees of freedom in order parameter
      integer isubd(7799)
* order parameter
* op(i,1,n) is the ith a component
* op(i,2,n) is the ith b component
* op(i,3,n) is the ith c component, etc.
      real op(8,8,7799)
* superspace group number for isotropy subgroup (1-775)
      integer ssgnumber(7799)
* superspace group number: 1.1, 2.1, 3.1, etc.
      character nstring(7799)*6
* symbol: P1(abg), P-1(abg), etc.
      character ssgstring(7799)*21
* basis vectors of lattice
* basis(i,j,n)=ith component of jth vector
      real basis(4,4,7799)
* origin of subgroup
      real origin(4,7799)

* open data file
      open(30,file='incommensurate_subgroups_data.txt')
* skip heading
      read(30,*)
      read(30,*)
      read(30,*)
* read each isotropy subgroup
      do m=1,7799
        read(30,*)n,ispacegroup(n),spacegroupstring(n),irnumber(n),
     $       irstring(n),kstring(n),kvectorstring(n),
     $       ((kvector(i,j,n),i=1,3),j=1,4),opstring(n),
     $       opvectorstring(n),irdim(n),isubd(n),
     $       ((op(i,j,n),i=1,irdim(n)),j=1,isubd(n)),
     $       ssgnumber(n),nstring(n),ssgstring(n),
     $       ((basis(i,j,n),i=1,4),j=1,4),(origin(i,n),i=1,4)
* fix fractions that have been rounded off
* k vector
        do j=1,4
          do i=1,3
            ng=nint(kvector(i,j,n)*12)
            if(abs(kvector(i,j,n)*12-ng).gt.0.01)then
              write(6,*)'kvec',kvector(i,j,n)
              stop
            endif
            kvector(i,j,n)=1.0*ng/12
          enddo
        enddo
* order parameter
        do j=1,isubd(n)
          do i=1,irdim(n)
            if(abs(op(i,j,n)).lt.0.01)then
              op(i,j,n)=0
            else
              sign=op(i,j,n)/abs(op(i,j,n))
              ng=nint(op(i,j,n)**2*12)
              if(abs(op(i,j,n)**2*12-ng).gt.0.01)then
                write(6,*)'op',op(i,j,n)
                stop
              endif
              op(i,j,n)=sign*sqrt(1.0*ng/12)
            endif
          enddo
        enddo
* basis
        do j=1,4
          do i=1,4
            ng=nint(basis(i,j,n)*12)
            if(abs(basis(i,j,n)*12-ng).gt.0.01)then
              write(6,*)'basis',basis(i,j,n)
              stop
            endif
            basis(i,j,n)=1.0*ng/12
          enddo
        enddo
* origin
        do i=1,4
          ng=nint(origin(i,n)*48)
          if(abs(origin(i,n)*48-ng).gt.0.01)then
            write(6,*)'origin',origin(i,n)
            stop
          endif
          origin(i,n)=1.0*ng/48
        enddo
      enddo
      end