anchor.diff:                          if (all(lvecs(:,ivec) == 0)) then
                               ! Nearest neighbours within unit cell.
                               call set_orb(tmat(:,i),j)
  -                            call set_orb(tmat(:,i+1),j+1)
  +                            if (isystem == 2) call set_orb(tmat(:,i+),j+1)
                           else if (.not. finite_cluster) then ! if we want inf. lattice
                               ! Nearest neighbours due to periodic boundaries.
                               call set_orb(tmat(:,j),i)
  -                            call set_orb(tmat(:,j+1),i+1) 
  +                            if (isystem == 2) call set_orb(tmat(:,j+1),i+1) 
                               ! else we just want connections to other cells to
                               ! stay as 0 
                           end if        
  
positive.code: integer, allocatable :: connected_sites(:,:) ! (2ndim, nbasis)
  
  ! True if any site is its own periodic image.
  ! This is the case if one dimension (or more) has only one site per crystal
  ! cell.  If so then the an orbital can incur a kinetic interaction with itself.
  ! This is the only way that the integral < i | T | i >, where i is a basis
  ! function centred on a lattice site, can be non-zero.
  
negative.code: 
          tmat = 0
          connected_orbs = 0
          
          ! For the Hubbard model, each orbital can have spin up or down, so 
          ! basis_fns(i) refers to alternating alpha and beta orbitals.
          ! In the do loop we therefore loop over every *second* orbital (because  
          ! spin must be the same for orbitals to be connected in this case).
          ! For Heisenberg, we just want to loop over every component of
          ! basis_fns, so we set isystem = 1
          if (system_type == heisenberg) then
              isystem = 1
          else
              isystem = 2
          endif
              
  
          ! Form all lattice vectors
          select case(ndim)
  