init_sites_triangular Subroutine

    subroutine init_sites_triangular(this)
        class(triangular) :: this
        integer :: temp_array(this%length(1), this%length(2))
        integer :: down(this%length(1), this%length(2))
        integer :: up(this%length(1), this%length(2))
        integer :: left(this%length(1), this%length(2))
        integer :: right(this%length(1), this%length(2))
        integer :: lu(this%length(1), this%length(2))
        integer :: rd(this%length(1), this%length(2))
        integer :: i, temp_neigh(6), x, y
        integer, allocatable :: neigh(:), order(:)
        character(*), parameter :: this_routine = "init_sites_triangular"

        ASSERT(this%get_nsites() >= 4)

        allocate(order(this%get_nsites()), source = 0)
        if (t_input_order .and. this%t_bipartite_order) then
            order = orbital_order
        else
            order = [(i, i = 1, this%get_nsites())]
        end if

        temp_array = reshape([(order(i), i=1, this%get_nsites())], this%length)

        up = cshift(temp_array, -1, 1)
        down = cshift(temp_array, 1, 1)
        right = cshift(temp_array, 1, 2)
        left = cshift(temp_array, -1, 2)
        lu = cshift(up, -1, 2)
        rd = cshift(down, 1, 2)

        if (this%is_periodic()) then
            do i = 1, this%get_nsites()

                x = mod(i - 1, this%length(1)) + 1
                y = (i - 1) / this%length(1) + 1

                temp_neigh = [up(x, y), down(x, y), left(x, y), right(x, y), lu(x, y), rd(x, y)]

                neigh = sort_unique(temp_neigh)

                this%sites(order(i)) = site(order(i), size(neigh), neigh)

                deallocate(neigh)

            end do
        else
            call stop_all(this_routine, &
                          "closed boundary conditions not yet implemented for triangular lattice!")
        end if

    end subroutine init_sites_triangular