function combine_spin_basis(n_orbs, n_first, n_second, n_total, first_basis, &
t_sort, n_doubles) result(spin_basis)
! function to combine a already given spin-basis with the second one
! with the additional option to specify the number of doubly occupied
! determinants..
! n_orbs ... is the number of spatial orbitals
! n_first ... number of spins, with which first_basis was created
! n_second ... number of spins with which basis is combined now
! n_total ... total size of basis for given number of doubly occ. sites!
! first_basis . already contructed one-spin basis
! t_sort ... should the basis be sorted. (after combination i am not sure it will be sorted..
! n_doubles .. (optional:) number of doubly occupied sites(has to fit with n_total!)
integer, intent(in) :: n_orbs, n_first, n_second, n_total
integer(n_int), intent(in) :: first_basis(:)
logical, intent(in) :: t_sort
integer, intent(in), optional :: n_doubles
integer(n_int) :: spin_basis(0:0, n_total)
character(*), parameter :: this_routine = "combine_spin_basis"
#ifdef DEBUG_
integer, allocatable :: n_dets(:)
#endif
integer :: n_doub, i, j, n_remain, n
integer(n_int), allocatable :: second_basis(:)
! the half-filled case is the most easy one.. should we treat it
! special? maybe..
! ASSERT(n_first >= n_second)
#ifdef DEBUG_
! be sure that the provided n_total fits with the n_doubles if
! provided
n_dets = calc_n_double(n_orbs, n_first, n_second)
if (present(n_doubles)) then
ASSERT(n_dets(n_doubles + 1) == n_total)
else
! n_doubles = 0 is the default
ASSERT(n_dets(1) == n_total)
end if
#endif
! default value of n_doubles = 0
if (present(n_doubles)) then
n_doub = n_doubles
else
n_doub = 0
end if
! do the n_double = 0 and half-filled case first, thats the
! easiest and only necessary for now..
select case (n_doub)
case (0)
if (n_first + n_second == n_orbs) then
! if (n_first == n_orbs) then
!
! else
! half-filled case: easiest
! since ms is not important, just convert the already
! calculated spin-basis to a spin-orbital basis:
! 0011 -> 00 00 01 01 eg.
do i = 1, n_total
! do have it ordered -> first set the beta spins on the left
spin_basis(:,i) = set_alpha_beta_spins(first_basis(i), n_orbs, .false.)
! and combine it with the alpha spins..
spin_basis(:,i) = ieor(spin_basis(:,i), set_alpha_beta_spins(not(first_basis(i)), n_orbs, .true.))
end do
! end if
else
! we have to distribute the n_second remaining spins
! across the n_orbs - n_first orbitals, so there are:
n_remain = int(choose_i64(n_orbs - n_first, n_second))
! states per first_basis states
ASSERT(n_remain * size(first_basis) == n_total)
second_basis = create_one_spin_basis(n_orbs - n_first, n_second)
! i want to do some sort of tensor product here:
! |0011> x |01> = |00 10 01 01>
! |0011> x |10> = |10 00 01 01>
n = 1
do i = 1, size(first_basis)
do j = 1, size(second_basis)
spin_basis(:,n) = open_shell_product(first_basis(i), second_basis(j), n_orbs)
n = n + 1
end do
end do
end if
case default
! i should not need a select case here.. this gets to
! complicated..
! it is similar to non-half filled actually..
if (n_first + n_second == n_orbs) then
! then we have the half-filled case..
else
end if
! not yet implemented
call stop_all(this_routine, "not yet implemented!")
end select
! do i want to sort it??
if (t_sort) then
call sort(spin_basis)
end if
end function combine_spin_basis