subroutine init_sites_cube(this)
class(cube) :: this
character(*), parameter :: this_routine = "init_sites_cube"
integer :: temp_array(this%length(1), this%length(2), this%length(3))
integer :: i, x, y, z, temp_neigh(6)
integer :: up(this%length(1), this%length(2), this%length(3))
integer :: down(this%length(1), this%length(2), this%length(3))
integer :: left(this%length(1), this%length(2), this%length(3))
integer :: right(this%length(1), this%length(2), this%length(3))
integer :: in(this%length(1), this%length(2), this%length(3))
integer :: out(this%length(1), this%length(2), this%length(3))
integer, allocatable :: neigh(:)
! minumum cube size:
ASSERT(this%get_nsites() >= 8)
! encode the lattice with the fortran intrinsic ordering of matrices
temp_array = reshape([(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)
in = cshift(temp_array, -1, 3)
out = cshift(temp_array, 1, 3)
if (this%is_periodic()) then
do i = 1, this%get_nsites()
! find conversion to 3D matrix indices..
! i am not yet sure about that: but seems to work
x = mod(i - 1, this%length(1)) + 1
z = (i - 1) / (this%length(1) * this%length(2)) + 1
y = mod((i - 1) / this%length(1), this%length(2)) + 1
temp_neigh = [up(x, y, z), down(x, y, z), left(x, y, z), right(x, y, z), &
in(x, y, z), out(x, y, z)]
neigh = sort_unique(temp_neigh)
this%sites(i) = site(i, size(neigh), neigh)
deallocate(neigh)
end do
else
call stop_all(this_routine, &
"closed boundary conditions not yet implemented for cubic lattice!")
end if
end subroutine init_sites_cube