subroutine apply_perturbation_array(perturbs, ndets, dets_in, dets_out, phase)
use bit_rep_data, only: NIfTot
use bit_reps, only: add_ilut_lists
use FciMCData, only: MaxWalkersPart
type(perturbation), intent(in) :: perturbs(:)
integer, intent(inout) :: ndets
integer(n_int), intent(in) :: dets_in(0:, :) ! First dimension must be 0:NIfTot
integer(n_int), intent(out) :: dets_out(0:, :) ! First dimension must be 0:NIfTot
real(dp), intent(in), optional :: phase(:) ! Phase factors of the perturbation operators
! size of phase must be equal to that of perturbs
integer :: i, ndets_init, ndets_pert_1, ndets_pert_2, ierr
integer(n_int), allocatable :: temp_dets_1(:, :), temp_dets_2(:, :)
ndets_init = ndets
if (size(perturbs) > 1) then
! Just allocate these arrays to be the same size as the CurrentDets
! array.
allocate(temp_dets_1(0:NIfTot, MaxWalkersPart), stat=ierr)
allocate(temp_dets_2(0:NIfTot, MaxWalkersPart), stat=ierr)
! Apply the first perturbation.
ndets_pert_1 = ndets_init
! Note that ndets_pert_1 is altered by this routine.
if (present(phase)) then
call apply_perturbation(perturbs(1), ndets_pert_1, dets_in, temp_dets_1, phase(1))
else
call apply_perturbation(perturbs(1), ndets_pert_1, dets_in, temp_dets_1)
end if
do i = 2, size(perturbs)
ndets_pert_2 = ndets_init
! Note that ndets_pert_2 is altered by this routine.
if (present(phase)) then
call apply_perturbation(perturbs(i), ndets_pert_2, dets_in, temp_dets_2, phase(i))
else
call apply_perturbation(perturbs(i), ndets_pert_2, dets_in, temp_dets_2)
end if
call add_ilut_lists(ndets_pert_1, ndets_pert_2, .false., temp_dets_1, temp_dets_2, dets_out, ndets)
! If we still have more perturbations to apply, copy the result
! to temp_dets_1. Else, exit the final result in dets_out.
if (i /= size(perturbs)) then
ndets_pert_1 = ndets
! this is very inefficient, but it is only performed at the beginning
! so leave it for now
temp_dets_1(:, 1:ndets) = dets_out(:, 1:ndets)
end if
end do
deallocate(temp_dets_1)
deallocate(temp_dets_2)
else if (size(perturbs) == 1) then
! Simply apply the single perturbation using the final input and
! output arrays.
! Ignore any given phase argument
call apply_perturbation(perturbs(1), ndets, dets_in, dets_out)
end if
end subroutine apply_perturbation_array