anchor.diff:  real(dp) :: Ainv(3,3), T(3,3), gtemp(3)
   logical err, mapped
   
  -open(18,file="debug_find_labeling.out")
  -write(18,'(3/,"Takes a list of vectors, maps them into the group, determines the labeling",2/)')
  +!open(18,file="debug_find_labeling.out")
  +!write(18,'(3/,"Takes a list of vectors, maps them into the group, determines the labeling",2/)')
   
   call matrix_inverse(A,Ainv,err)
   if(err) stop "Coplanar vectors in find_labeling_from_atom_basis"
  
positive.code: nAt = size(aBas,2)
  nD = size(dset,2) ! Number of d-vectors in the dset
  
  write(18,'("Diagonal entries of HNF: ",3(i3,1x))') diag
  write(18,'("Left transform:",/,3(3i3,1x,/),/)') transpose(L)
  write(18,'("parent lattice vectors (columns):",/,3(3(f7.3,1x),/),/)') transpose(A)
  write(18,'("parLatt inverse (columns):",/,3(3(f7.3,1x),/),/)') transpose(Ainv)
  do iAt = 1, nAt
     write(18,'("Atom #: ",i3,"   position: ",3(f7.3,1x))') iAt,aBas(:,iAt)
  enddo
  
  g = 0
  do iAt = 1, nAt ! Map each real space vector (atom position) into the group
  
negative.code:    p_with_d(1:3,pLindx:pLindx+n-1) = p
     p_with_d(4,  pLindx:pLindx+n-1) = iD
  enddo
  !do i = 1, 4
  !   write(17,'(400(i2,1x))') p_with_d(i,:)
  !enddo
  
  call matrix_inverse(matmul(L,Ainv),T)
  !write(18,'("Map from group:",/,3(3(f7.3,1x),/),/)') transpose(T)
  
  ! If we want to do this, we need the supercell vectors
  !do iAt = 1, nAt
  