anchor.diff:  
           tmat = 0
           connected_orbs = 0
  +
  +        ! Form all lattice vectors
  +        select case(ndim)
  +        case(1)
  +            do i = -1, 1
  +                lvecs(:,i+2) = i*lattice(:,1)
  +            end do
  +        case(2)
  +            do i = -1, 1
  +                do j = -1, 1
  +                    lvecs(:,j+2+3*(i+1)) = i*lattice(:,1) + j*lattice(:,2)
  +                end do
  +            end do
  +        case(3)
  +            do i = -1, 1
  +                do j = -1, 1
  +                    do k = -1, 1
  +                        lvecs(:,k+2+3*(j+1)+9*(i+1)) = i*lattice(:,1) + j*lattice(:,2) + k*lattice(:,3)
  +                    end do
  +                end do
  +            end do
  +        end select
   
           ! Construct how the lattice is connected.
           do i = 1, nbasis-1, 2
  
positive.code: 
          use basis, only: nbasis, bit_lookup, basis_length, basis_fns, set_orb
          use system, only: lattice, ndim, box_length
  
          integer :: i, j, ierr, pos, ind, ivec
          integer :: r(ndim)
  
          t_self_images = any(abs(box_length-1.0_p) < depsilon)
  
  
negative.code:                 ! filled in automatically.
                  ! All matrix elements between different spins are zero
                  ! Allow j=i in case i is its own periodic image.
                  r = basis_fns(i)%l - basis_fns(j)%l
                  do ivec = 1, 3**ndim
                      if (sum(abs(r-lvecs(:,ivec))) == 1) then
                          ! i and j are on sites which are nearest neighbours 
                          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)
                          else
                              ! Nearest neighbours due to periodic boundaries.
                              call set_orb(tmat(:,j),i)
                              call set_orb(tmat(:,j+1),i+1)
                          end if
                          call set_orb(connected_orbs(:,i),j)
                          call set_orb(connected_orbs(:,j),i)
                          call set_orb(connected_orbs(:,i+1),j+1)
  