init_sites_hexagonal Subroutine

    subroutine init_sites_hexagonal(this)
        ! the way to create the neighboring indices is based on the convention
        ! that the lattice is interpreted in such a way: the unit cell is:
        !  __
        ! /  \   with encoding: 1 5
        ! \__/                  2 6
        ! /  \                  3 7
        !                       4 8

        ! which leads to an lattice:
        ! __/  \__/
        !   \__/  \_
        ! _ /  \__/
        !   \__/  \_
        ! __/  \__/
        !   \__/  \_

        ! so the up and down neighbors are easy to do again.
        ! just the left and right are different: now each site only either
        ! has left or right, not both and it alternates
        ! but this also depends on the colum of encoding one is in:
        !   1 - 5    9 - 13
        ! - 2   6 - 10   14 -
        !   3 - 7   11 - 15
        ! - 4   8 - 12   16 -
        class(hexagonal) :: this
        integer :: temp_array(4 * this%length(1), 2 * this%length(2))
        integer :: up(4 * this%length(1), 2 * this%length(2))
        integer :: down(4 * this%length(1), 2 * this%length(2))
        integer :: right(4 * this%length(1), 2 * this%length(2))
        integer :: left(4 * this%length(1), 2 * this%length(2))

        integer :: i, temp_neigh(3), x, y, special
        integer, allocatable :: neigh(:), order(:)
        character(*), parameter :: this_routine = "init_sites_hexagonal"

        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())], &
                             [4 * this%length(1), 2 * this%length(2)])

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

        if (this%is_periodic()) then
            do i = 1, this%get_nsites()
                ! columns and rows:
                x = mod(i - 1, 4 * this%length(1)) + 1
                y = (i - 1) / (4 * this%length(1)) + 1

                if (is_odd(y)) then
                    if (is_odd(i)) then
                        ! every odd number in a odd column has a right neighbor
                        special = right(x, y)
                    else
                        ! otherwise left
                        special = left(x, y)
                    end if
                else
                    ! for even columns it is the other way around
                    if (is_odd(i)) then
                        special = left(x, y)
                    else
                        special = right(x, y)
                    end if
                end if
                temp_neigh = [up(x, y), down(x, y), special]

                neigh = sort_unique(temp_neigh)

                this%sites(order(i)) = site(order(i), size(neigh), neigh)
            end do
        else
            call stop_all(this_routine, &
                          "closed boundary conditions not yet implemented for hexagonal lattice!")

        end if

    end subroutine init_sites_hexagonal