anchor.diff:  !integer diag(3)
   ! debug
   
  -
   allocate(label(1,size(dset,2)))
   allocate(digit(size(dset,2)))
   label = 1
   digit = 1
   
  -
  -open(17,file="debug_gspace_rep.out")
  +! open(17,file="debug_gspace_rep.out")
   nAt = size(aBas,2)
  -write(17,'("Number of atoms: ",i3)') nAt
  +! write(17,'("Number of atoms: ",i3)') nAt
   if(nAt/=size(aTyp)) stop "Input to get_gspace_representation is inconsistent: nAt"
   
   
  
positive.code: call find_equivalent_labelings(labeling,LattRotList,pLabel)
  
  nuq = size(pLabel,1)
  write(17,'(/,"Number of unique labelings: ",i5)') nuq
  do iuq = 1, nuq
     write(17,'("uq Labeling # :",i3,5x,"labeling:",1x,200(i1,1x))') iuq,pLabel(iuq,:)
  enddo
  close(17)
  END SUBROUTINE get_gspace_representation
  
  !***************************************************************************************************
  
negative.code:   ! rotationally-equivalent. So the input and output lists are only one element
    ! long (in the last index). E.g., HNFin is a 3x3x1 array of integers (original lattice HNF)
  call remove_duplicate_lattices(HNFin,LatDim,pLV,dset,dRotList,HNFout,fixOp,LattRotList,sLVlist,label,digit,eps)
  ! write(17,'("Number of symmetry operations that fix the superlattice: ",i3,/)') size(fixOp(1)%rot,3)
  ! write(17,'(200(3(3f7.3,1x,/),"shift:",3(f7.3,1x),//))') &
  !     ((fixOp(1)%rot(j,:,iOp),j=1,3),fixOp(1)%shift(:,iOp),iOp=1,size(fixOp(1)%rot,3))
    ! This routine gets the SNF form of the HNF, returns a list of permutations effected by the
    ! rotation symmetries,
  call get_SNF(HNFout,L,SNFlist,R,LattRotList,uqSNF,SNFlabel,fixOp)
  