subroutine init_sites_rect(this)
class(rectangle) :: this
character(*), parameter :: this_routine = "init_sites_rect"
integer :: i, temp_neigh(4), x, y
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, allocatable :: neigh(:)
integer :: k_vec(3), r_vec(3)
integer, allocatable :: order(:)
! this is the important routine..
! store lattice like that:
! 1 4 7
! 2 5 8
! 3 6 9
ASSERT(this%get_nsites() >= 4)
! use cshift intrinsic of fortran..
! how do i efficiently set that up?
if (this%t_bipartite_order) then
if (this%length(1) /= this%length(2)) then
if (this%length(2) /= 2) then
call stop_all(this_routine, &
"ladder bipartite ordering is implemented with Ly == 2)")
end if
allocate(order(this%get_nsites()), source = 0)
if (t_input_order) then
order = orbital_order
else
do i = 1, this%length(1)/2
order(2*i-1) = i
order(2*i) = this%length(1) + i
end do
do i = this%length(1)/2 + 1, this%length(1)
order(2*i-1) = this%length(1) + i
order(2*i) = i
end do
end if
else
if (this%get_nsites() == 16) then
allocate(order(16), source = 0)
if (t_input_order) then
order = orbital_order
else
order = [ 1, 9, 2, 10, &
11, 3, 12, 4, &
5, 13, 6, 14, &
15, 7, 16, 8]
end if
else if (this%get_nsites() == 36) then
allocate(order(36), source = 0)
if (t_input_order) then
order = orbital_order
else
order = [ 1, 19, 2, 20, 3, 21, &
22, 4, 23, 5, 24, 6, &
7, 25, 8, 26, 9, 27, &
28, 10, 29, 11, 30, 12, &
13, 31, 14, 32, 15, 33, &
34, 16, 35, 17, 36, 18]
end if
else
if (t_input_order) then
order = orbital_order
else
call stop_all(this_routine, &
"bipartite order for square only implemented for 4x4! and 6x6 for now!")
end if
end if
end if
else
allocate(order(this%get_nsites()), source = [(i, i = 1, this%get_nsites())])
end if
temp_array = reshape( order, 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)
if (this%is_periodic()) then
do i = 1, this%get_nsites()
! create the neighbor list
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)]
neigh = sort_unique(temp_neigh)
k_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
r_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
this%sites(order(i)) = site(order(i), size(neigh), neigh, k_vec, r_vec)
deallocate(neigh)
end do
else if (this%is_periodic(1)) then
! only periodic in x-direction
do i = 1, this%get_nsites()
x = mod(i - 1, this%length(1)) + 1
y = (i - 1) / this%length(1) + 1
! now definetly always take the left and right neighbors
! but up and down if we are not jumping boundaries
if (x == 1) then
! dont take upper neighbor -> just repeat an neighbor,
! which will get removed from unique
temp_neigh = [down(x, y), down(x, y), left(x, y), right(x, y)]
else if (x == this%length(1)) then
temp_neigh = [up(x, y), up(x, y), left(x, y), right(x, y)]
else
! take all
temp_neigh = [up(x, y), down(x, y), left(x, y), right(x, y)]
end if
neigh = sort_unique(temp_neigh)
k_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
r_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
this%sites(order(i)) = site(order(i), size(neigh), neigh, k_vec, r_vec)
deallocate(neigh)
end do
else if (this%is_periodic(2)) then
! only periodic in the y-direction
do i = 1, this%get_nsites()
x = mod(i - 1, this%length(1)) + 1
y = (i - 1) / this%length(1) + 1
! now definetly always take the left and right neighbors
! but up and down if we are not jumping boundaries
if (y == 1) then
! dont take upper neighbor -> just repeat an neighbor,
! which will get removed from unique
temp_neigh = [up(x, y), down(x, y), right(x, y), right(x, y)]
else if (y == this%length(2)) then
temp_neigh = [up(x, y), down(x, y), left(x, y), left(x, y)]
else
! take all
temp_neigh = [up(x, y), down(x, y), left(x, y), right(x, y)]
end if
neigh = sort_unique(temp_neigh)
k_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
r_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
this%sites(order(i)) = site(order(i), size(neigh), neigh, k_vec, r_vec)
deallocate(neigh)
end do
else
! non-periodic
do i = 1, this%get_nsites()
x = mod(i - 1, this%length(1)) + 1
y = (i - 1) / this%length(1) + 1
! now definetly always take the left and right neighbors
! but up and down if we are not jumping boundaries
if (x == 1 .and. y == 1) then
! dont take upper neighbor -> just repeat an neighbor,
! which will get removed from unique
temp_neigh = [down(x, y), down(x, y), right(x, y), right(x, y)]
else if (y == this%length(2) .and. x == this%length(1)) then
temp_neigh = [up(x, y), up(x, y), left(x, y), left(x, y)]
else if (x == this%length(1) .and. y == 1) then
temp_neigh = [up(x, y), up(x, y), right(x, y), right(x, y)]
else if (y == this%length(2) .and. x == 1) then
temp_neigh = [down(x, y), down(x, y), left(x, y), left(x, y)]
else if (x == 1) then
temp_neigh = [left(x, y), right(x, y), down(x, y), down(x, y)]
else if (x == this%length(1)) then
temp_neigh = [left(x, y), right(x, y), up(x, y), up(x, y)]
else if (y == 1) then
temp_neigh = [right(x, y), right(x, y), up(x, y), down(x, y)]
else if (y == this%length(2)) then
temp_neigh = [left(x, y), left(x, y), up(x, y), down(x, y)]
else
! take all
temp_neigh = [up(x, y), down(x, y), left(x, y), right(x, y)]
end if
neigh = sort_unique(temp_neigh)
k_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
r_vec = [x - (this%length(1) + 1) / 2, y - (this%length(2) + 1) / 2, 0]
this%sites(order(i)) = site(order(i), size(neigh), neigh, k_vec, r_vec)
deallocate(neigh)
end do
end if
end subroutine init_sites_rect