function create_neel_state(ilut_neel) result(neel_state)
! probably a good idea to have a routine which creates a neel state
! (or one of them if it is ambigous)
integer(n_int), intent(out), optional :: ilut_neel(0:NIfTot)
integer :: neel_state(nel)
character(*), parameter :: this_routine = "create_neel_state"
integer(n_int), allocatable :: tmp_ilut(:)
integer :: i, j, k, l, spin, ind
! make this independent of the lattice mod, to call it as early
! as possible in the initialization
! also assert that we have at most half-filling.. otherwise it does
! not make so much sense to do this
! there is no neel state for all lattice, but atleast something
! similar to it..
! atleast for the lattice where there is a neel state, i should
! create it automaticall and for the other a similar one
! the lattice should already be intialized
! i was thinking of putting this functionality into the
! lattice_mod, but i guess i should keep it seperate, since it
! actually has to do more with the system or?
neel_state = 0
select case (lattice_type)
case ('chain')
neel_state = create_neel_state_chain()
case ('square', 'rectangle', 'triangle', 'triangular', 'hexagonal', 'kagome')
! check if length_x is mod 2
if (mod(length_x, 2) == 0) then
neel_state = create_neel_state_chain()
! and flip the spins in every other column
do i = 1, length_y, 2
! loop over every second column
if (i * length_x >= nel) exit
do j = 1, length_x
! whats the total index?
ind = i * length_x + j
if (ind > nel) exit
! flip the spins.. this should be way easier..
if (is_beta(neel_state(ind))) then
neel_state(ind) = neel_state(ind) + 1
else
neel_state(ind) = neel_state(ind) - 1
end if
end do
end do
else
! here it is easy it is just like the chain case
neel_state = create_neel_state_chain()
end if
case ('cube')
! not yet implemented
ASSERT(.false.)
case ('tilted', 'tilted-square', 'square-tilted')
! do is similar to the actual construction of the lattice
! but way nicer as this stuff below..
k = 0
l = 1
spin = 1
do i = -length_x + 1, 0
do j = -k, k
neel_state(l) = spin
l = l + 1
if (l > nel) exit
if (is_beta(spin)) then
spin = spin + 3
else
spin = spin + 1
end if
end do
k = k + 1
if (l > nel) exit
! in the first half we need the same starting spin
spin = spin - 1
end do
k = k - 1
! finished already?
if (l > nel) return
spin = spin + 1
do i = 1, length_x
do j = -k, k
neel_state(l) = spin
l = l + 1
if (l > nel) exit
if (is_beta(spin)) then
spin = spin + 3
else
spin = spin + 1
end if
end do
k = k - 1
if (l > nel) exit
spin = spin + 1
end do
case default
call stop_all(this_routine, "unknown lattice type!")
end select
if (tHPHF) then
! i have to get the relevant HPHF determinant
nbasis = 2 * lat%get_nsites()
call init_bit_rep()
allocate(tmp_ilut(0:niftot), source = 0_n_int)
call EncodeBitDet(neel_state, tmp_ilut)
tmp_ilut = return_hphf_sym_det(tmp_ilut)
call decode_bit_det(neel_state, tmp_ilut)
deallocate(tmp_ilut)
end if
if (present(ilut_neel)) then
call EncodeBitDet(neel_state, ilut_neel)
end if
end function create_neel_state