anchor.diff:  allocate(aTypTemp(nAt),aBasTemp(3,nAt))
   aTypTemp = aTyp; aBasTemp = aBas
   pLVtemp = sLV 
  -call make_primitive(pLVtemp,aTypTemp,aBasTemp,.false.,eps)
  +call make_primitive(pLVtemp,aTypTemp,aBasTemp,.false.,eps,removed_=removedAtoms)
   if(nAt/=size(aTypTemp)) then; ! the structure was reduced by make_primitive
      ! write(17,'(/,"ERROR: The input structure wasn''t primitive")')
      ! write(17,'("atom type: ",20(i2,1x))') aTypTemp(:)
      ! write(17,'("number of atoms: ",i2,5x," size of aTyp:",i2)') nAt, size(aTypTemp)
  -   stop "ERROR: input structure for get_HNF_of_derivative_structure was not primitive";endif
  +   !stop "ERROR: input structure for get_HNF_of_derivative_structure was not primitive";endif
  +  write(*,'(/,A)') "WARNING: input structure for get_HNF_of_derivative_structure was not primitive"; 
  +  allocate(hnf(1,1,1)); HNF = 0
  +  return
  +endif
   
   ! Now make every atom the same and apply make_primitive to find the parent cell
  -aTypTemp = 1; pLVtemp = sLV
  +! (but don't change the labeling of fixed occupation sites)
  +aTypTemp = 1
  +aTypTemp(fixedOcc_) = aTyp(fixedOcc_)
  +pLVtemp = sLV
   call make_primitive(pLVtemp,aTypTemp,aBasTemp,.false.,eps)
   ! write(17,'("Parent lattice of input superlattice (columns): ",/,3(3(f8.4,1x),/))') (pLVtemp(iAt,:),iAt=1,3)
   
   call matrix_inverse(pLVtemp,parLattTest,err)
   if(err) stop "Problem inverting parent lattice basis"
  -if (nD/=size(aTypTemp))then;print*,"Number of parent lattice sites in " !//trim(adjustl(sfname))
  +if (nD/=size(aTypTemp))then;print*,"Number of parent lattice sites" !in " //trim(adjustl(sfname))
      print *,"isn't the same as the structure being checked";
      print *,"size(aTypTemp)",size(aTypTemp),"nD",nD;stop;endif
   ! If we pass this test, the two cell bases are equivalent even if not equal. From this point on, use
   ! the one from the struct_enum.out file. (that one is pLV, not pLVtemp)
   
   
  -allocate(dsetStruc(3,size(aTypTemp))); dsetStruc = aBasTemp
  +allocate(dsetStruc(3,nD)); dsetStruc = aBasTemp
   call matrix_inverse(pLV,pLVinv,err) 
   if(err) stop "matrix inverse failed in get_HNF_of_derivative_structure"
   ! write(17,'("d-set: ",/,200(3(f8.4,1x),/))') (dsetStruc(:,iAt),iAt=1,size(dsetStruc,2))
  
positive.code: !* Check that the site basis for both the enum file and POSCAR are consistent.
  !  They might not have the same origin so try all distinct origin shifts and see if there is any
  !  shift that maps all of the sites in one case to all of the sites in the other
  do iAt = 1, nD
     !write(*,'("iAt:",i3,"  EnumBas ",3(F8.4,1x))') iAt,EnumBas(:,iAt)
     !write(*,'("iAt:",i3,"  dset ",3(F8.4,1x))') iAt,dset(:,iAt)
     diff = dsetStruc(:,iAt)-dset(:,1)
     !write(*,'("iAt:",i3,"  diff ",3(F8.4,1x))') iAt,dset(:,iAt)-EnumBas(:,1)
     do iE = 1, nD
        mapped=.false.
        do iP = 1, nD
           testsite = dset(:,iP)-diff
           !write(*,'("iP:",i3,"  testsite",3(F8.4,1x))') iP,testsite
           call bring_into_cell(testsite,pLVinv,pLV,eps) ! Make sure we stay in the same cell        
           !write(*,'("iP:",i3,"  testsite",3(F8.4,1x))') iP,testsite
           if (equal(testsite,dset(:,iE),eps)) then ! this site maps to one in other lattice
              mapped=.true.;!print *,"mapped" 
              exit
           endif
        enddo
        if(.not. mapped) exit
     enddo
     ! If we get to here and mapped is true, then the iE loop concluded without mapped ever coming up
     ! false in the iP loop. If that is the case, then the sites in EnumBas all had coincident sites
     ! in aBas. That means that the parent lattices (POSCAR and enum) are equivalent even if they
     ! don't have the same origin.
     if (mapped) exit
  enddo
  if (.not. mapped) stop "The lattice sites of the two parent lattices are not coincident"
  if (.not. equal(diff,0._dp,eps)) then
     ! write(17,'("The atomic sites of the parent lattice of the input structure have been shifted by")')
     ! write(17,'(3(f8.4,1x),/)') diff
     dset = dset - spread(diff,2,size(aBas(:,2)))
     do iAt = 1, nD
        call bring_into_cell(dset(:,iAt),pLVinv,pLV,eps)
        ! write(17,'("new pBas: ",3(f8.4,1x))') dset(:,iAt)
  
negative.code: !integer diag(3)
  ! debug
  
  !write(*,'("d-set: ",/,200(3(f7.3,1x),/))') (dset(:,iAt),iAt=1,size(dset,2))
  
  allocate(label(1,size(dset,2)))
  allocate(digit(size(dset,2)))
  label = 1
  