anchor.diff:          call check_allocate('tmat',basis_length*nbasis,ierr)
           allocate(connected_orbs(basis_length,nbasis), stat=ierr)
           call check_allocate('connected_orbs',basis_length*nbasis,ierr)
  +        allocate(connected_sites(2*ndim,nbasis), stat=ierr)
  +        call check_allocate('connected_sites',basis_length*2*ndim,ierr)
   
           tmat = 0
           connected_orbs = 0
  
positive.code:         ! Initialise real space Hubbard model: find and store the matrix
          ! elements < i | T | j > where i and j are real space basis functions.
  
          use basis, only: nbasis, bit_lookup, basis_length, basis_fns, set_orb
          use system, only: lattice, ndim, box_length
          use bit_utils
          use checking, only: check_allocate
          use errors, only: stop_all
  
          integer :: i, j, k, ierr, pos, ind, ivec
          integer :: r(ndim)
  
          integer :: lvecs(ndim, 3**ndim)
  
negative.code:                             ! Nearest neighbours within unit cell.
                              call set_orb(tmat(:,i),j)
                              call set_orb(tmat(:,i+1),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) 
                              ! else we just want connections to other cells to
                              ! stay as 0 
                          end if        
                         
                          ! If we only want a discrete molecule and the lattice
                          ! vector connecting the 2 sites is the 0-vector then the
                          ! 2 sites are connected in a unit cell and thus are
                          ! actually connected. (If they "connect" accross cell
                          ! boundaries then they are not connected for a single
                          ! molecule).
                          if ( (finite_cluster .and. all(lvecs(:,ivec) == 0)) .or. &
                               .not. finite_cluster) then
                              if (i /= j) then
                                  ! connected_orbs does not contain self-connections 
                                  ! due to the periodic boundary conditions.
                                  call set_orb(connected_orbs(:,i),j)
                                  call set_orb(connected_orbs(:,i+1),j+1)                      
                                  call set_orb(connected_orbs(:,j),i)
                                  call set_orb(connected_orbs(:,j+1),i+1)
                              end if
                          end if
                      end if
                  end do
              end do
          end do
  
          ! Decode connected_orbs to store list of connections.
          connected_sites = 0
          do i = 1, nbasis
              v = 0
              do ind = 1, basis_length
                  do pos = 0, i0_end
                      if (btest(connected_orbs(ind,i), pos)) then
                          v = v + 1
                          connected_sites(v, i) = basis_lookup(pos, ind)
                      end if
                  end do
              end do
          end do
  
      end subroutine init_real_space_hub
  