subroutine perform_multiplication(ras, classes, ras_strings, ras_iluts, ras_excit, vec_in, vec_out)
! In: ras - Details of the RAS space being used.
! In: classes - Details of the RAS classes for the above RAS space. This array (and also the ras
! variable) should be created by the the initialise_ras_space routine in the ras module.
! In: ras_strings - All of the alpha (and beta) strings for the RAS space.
! In: ras_iluts - All of the alpha (and beta) iluts for the RAS space.
! In: ras_excit - For each RAS string, store the addresses of strings which can be excited to
! by a single excitation. This array, and the two above, are created by the
! create_direct_ci_arrays routine in this module.
! In: vec_in - The input vector to be multiplied, in block form. A standard vector can be bought
! to this form by using the transfer_to_block_form routine in this module.
! Inout: vec_out - The resulting output vector after multiplication.
type(ras_parameters), intent(in) :: ras
type(ras_class_data), intent(in) :: classes(ras%num_classes)
integer, intent(in) :: ras_strings(-1:tot_nelec, ras%num_strings)
integer(n_int), intent(in) :: ras_iluts(0:NIfD, ras%num_strings)
type(direct_ci_excit), intent(in) :: ras_excit(ras%num_strings)
type(ras_vector), intent(in) :: vec_in(ras%num_classes, ras%num_classes, 0:7)
type(ras_vector), intent(inout) :: vec_out(ras%num_classes, ras%num_classes, 0:7)
type(ras_factors) :: factors(ras%num_classes, 0:7)
real(dp), allocatable, dimension(:) :: alpha_beta_fac
integer, allocatable, dimension(:) :: r
real(dp), allocatable, dimension(:, :) :: c
integer, allocatable, dimension(:, :) :: excit_info
integer :: string_i(tot_nelec), string_j(tot_nelec)
integer(n_int) :: ilut_i(0:NIfD)
logical :: in_ras_space
real(dp) :: v
integer :: class_i, class_j, class_k, class_m
integer :: excited_class, excited_sym
integer :: par_1, par_2
integer :: i, j, k, l, m
integer :: excit_j, excit_k
integer :: ind_i, ind_j, ind_k, full_ind_j, full_ind_k
integer :: nras1, nras3, min_ind, max_ind
integer :: sym_i, sym_j, sym_k, sym_m
integer :: ex1(2), ex2(2)
integer :: BRR_ex1(2), BRR_ex2(2)
integer :: ierr
! Initialisation - allocate arrays.
do class_i = 1, ras%num_classes
do sym_i = 0, 7
allocate(factors(class_i, sym_i)%elements(1:classes(class_i)%num_sym(sym_i)), &
stat=ierr)
end do
end do
do class_i = 1, ras%num_classes
do j = 1, classes(class_i)%num_comb
class_j = classes(class_i)%allowed_combns(j)
do sym_i = 0, 7
sym_j = ieor(HFSym_ras, sym_i)
if (classes(class_i)%num_sym(sym_i) == 0 .or. &
classes(class_j)%num_sym(sym_j) == 0) cycle
vec_out(class_i, class_j, sym_i)%elements(:, :) = 0.0_dp
end do
end do
end do
allocate(alpha_beta_fac(ras%num_strings))
allocate(c(ras%num_strings, ras%num_strings))
allocate(r(ras%num_strings))
alpha_beta_fac = 0.0_dp
c = 0.0_dp
r = 0
allocate(excit_info(2, ras%num_strings))
excit_info = 0
! Beta and beta-beta (sigma_1) and alpha and alpha-alpha (sigma_2) contributions.
! See Eq. 20 for algorithm.
! Loop over all strings.
do i = 1, ras%num_strings
! Set factors array equal to zero.
call zero_factors_array(ras%num_classes, factors)
! Get info for string_i.
sym_i = ras_strings(-1, i)
class_i = ras_strings(0, i)
ind_i = i - ras%cum_classes(class_i) - classes(class_i)%cum_sym(sym_i)
! Loop over all single excitations from string_i (-> string_k).
do excit_k = 1, ras_excit(i)%nexcit
! The full address of string k.
full_ind_k = ras_excit(i)%excit_ind(excit_k)
sym_k = ras_strings(-1, full_ind_k)
class_k = ras_strings(0, full_ind_k)
par_1 = ras_excit(i)%par(excit_k)
ex1 = ras_excit(i)%orbs(:, excit_k)
BRR_ex1(1) = BRR(2 * ex1(1)) / 2
BRR_ex1(2) = BRR(2 * ex1(2)) / 2
! The shifted address (shifted so that the first string with this symmetry and
! in this RAS class has address 1).
ind_k = full_ind_k - ras%cum_classes(class_k) - classes(class_k)%cum_sym(sym_k)
! In Eq. 20, this term is sgn(kl)*h_{kl}.
factors(class_k, sym_k)%elements(ind_k) = factors(class_k, sym_k)%elements(ind_k) + &
par_1 * GetTMatEl(BRR_ex1(2) * 2, BRR_ex1(1) * 2)
! The following lines add in g_{kl} from Eq. 29.
do j = 1, ex1(2) - 1
factors(class_k, sym_k)%elements(ind_k) = factors(class_k, sym_k)%elements(ind_k) - &
par_1 * get_umat_el(BRR_ex1(2), BRR(2 * j) / 2, BRR(2 * j) / 2, BRR_ex1(1))
end do
if (ex1(2) > ex1(1)) then
factors(class_k, sym_k)%elements(ind_k) = factors(class_k, sym_k)%elements(ind_k) - &
par_1 * get_umat_el(BRR_ex1(2), BRR_ex1(2), BRR_ex1(2), BRR_ex1(1))
else if (ex1(2) == ex1(1)) then
factors(class_k, sym_k)%elements(ind_k) = factors(class_k, sym_k)%elements(ind_k) - &
0.5_dp * par_1 * get_umat_el(BRR_ex1(2), BRR_ex1(2), BRR_ex1(2), BRR_ex1(1))
end if
! Loop over all single excitations from string_k (-> string_j).
do excit_j = 1, ras_excit(full_ind_k)%nexcit
ex2 = ras_excit(full_ind_k)%orbs(:, excit_j)
! Only need to consider excitations where (ij) >= (kl) (see Eq. 28).
if ((ex2(2) - 1) * tot_norbs + ex2(1) < (ex1(2) - 1) * tot_norbs + ex1(1)) cycle
BRR_ex2(1) = BRR(2 * ex2(1)) / 2
BRR_ex2(2) = BRR(2 * ex2(2)) / 2
! The full address for string j.
full_ind_j = ras_excit(full_ind_k)%excit_ind(excit_j)
sym_j = ras_strings(-1, full_ind_j)
class_j = ras_strings(0, full_ind_j)
par_2 = ras_excit(full_ind_k)%par(excit_j)
! The shifted address (shifted so that the first string with this symmetry and
! in this RAS class has address 1).
ind_j = full_ind_j - ras%cum_classes(class_j) - classes(class_j)%cum_sym(sym_j)
! Avoid overcounting for the case that the indices are the same.
if (ex1(1) == ex2(1) .and. ex1(2) == ex2(2)) then
factors(class_j, sym_j)%elements(ind_j) = factors(class_j, sym_j)%elements(ind_j) + &
0.5_dp * par_1 * par_2 * get_umat_el(BRR_ex2(2), BRR_ex1(2), BRR_ex2(1), BRR_ex1(1))
else
factors(class_j, sym_j)%elements(ind_j) = factors(class_j, sym_j)%elements(ind_j) + &
par_1 * par_2 * get_umat_el(BRR_ex2(2), BRR_ex1(2), BRR_ex2(1), BRR_ex1(1))
end if
end do
end do
! The factors array has now been fully generated for string_i. Now we just have to
! add in the contribution to vec_out from this string.
! This performs the final line in Eq. 20.
! Loop over all classes connected to the class of string_i.
do j = 1, classes(class_i)%num_comb
class_j = classes(class_i)%allowed_combns(j)
! The *required* symmetry,
sym_j = ieor(HFSym_ras, sym_i)
! If there are no states in this class with the required symmetry.
if (classes(class_j)%num_sym(sym_j) == 0) cycle
! Loop over all classes connected to string_j.
do k = 1, classes(class_j)%num_comb
class_k = classes(class_j)%allowed_combns(k)
! The *required* symmetry, sym_k = ieor(HFSym_ras, sym_j) = sym_i.
sym_k = sym_i
! If there are no states in this class with the required symmetry.
if (classes(class_k)%num_sym(sym_k) == 0) cycle
! Add in sigma_2.
vec_out(class_j, class_i, sym_j)%elements(:, ind_i) = &
vec_out(class_j, class_i, sym_j)%elements(:, ind_i) + &
matmul(vec_in(class_j, class_k, sym_j)%elements(:, :), factors(class_k, sym_k)%elements)
! Add in sigma_1.
vec_out(class_i, class_j, sym_i)%elements(ind_i, :) = &
vec_out(class_i, class_j, sym_i)%elements(ind_i, :) + &
matmul(factors(class_k, sym_k)%elements, vec_in(class_k, class_j, sym_k)%elements(:, :))
end do ! Over all classes connected to string_j.
end do ! Over all classes connected to string_i.
end do ! Over all strings.
! Next, calculate the contribution from the alpha-beta term (sigma_3) (Eq. 23).
! Loop over all combinations of two spatial indices, (kl).
do k = 1, tot_norbs
do l = 1, tot_norbs
! Zero arrays.
r = 0
! The array c here is C' in Eq. 23.
c = 0.0_dp
ex1(1) = l
ex1(2) = k
! Store these, for quick access later.
BRR_ex1(1) = BRR(2 * l) / 2
BRR_ex1(2) = BRR(2 * k) / 2
! Loop over all strings (equivalent to J_{beta} in Eq. 23).
do i = 1, ras%num_strings
! Get info for string_i.
sym_i = ras_strings(-1, i)
class_i = ras_strings(0, i)
nras1 = classes(class_i)%nelec_1
nras3 = classes(class_i)%nelec_3
string_i = ras_strings(1:tot_nelec, i)
ilut_i = ras_iluts(:, i)
! This is the condition for (kl) to be an allowed excitation from the current string.
if (IsOcc(ilut_i, l) .and. &
(.not. (IsOcc(ilut_i, k) .and. k /= l))) then
! Temporarily set these for get_excit_details to use.
sym_j = sym_i
class_j = class_i
string_j = string_i
call get_excit_details(ex1, ras, nras1, nras3, string_j, sym_j, class_j, in_ras_space)
! If the excitation is to outside the RAS space, then we don't need to consider it.
if (in_ras_space) then
! Store the class and symmetry of the excited string for later, to save some speed.
excit_info(1, i) = class_j
excit_info(2, i) = sym_j
! The full address for string j.
ind_j = classes(class_j)%address_map(get_address(classes(class_j), string_j))
! The shifted address, so that the first string in this class will have address 1.
ind_j = ind_j - classes(class_j)%cum_sym(sym_j)
r(i) = get_single_parity(ilut_i, l, k)
! Construct array C' in Eq 23. To do so, loop over all connected strings.
do m = 1, classes(class_j)%num_comb
class_m = classes(class_j)%allowed_combns(m)
sym_m = ieor(HFSym_ras, sym_j)
if (classes(class_m)%num_sym(sym_m) /= 0) then
min_ind = ras%cum_classes(class_m) + classes(class_m)%cum_sym(sym_m) + 1
max_ind = min_ind + classes(class_m)%num_sym(sym_m) - 1
c(min_ind:max_ind, i) = real(r(i), dp) * vec_in(class_j, class_m, sym_j)%elements(ind_j, :)
end if
end do
end if
end if
end do
! Loop over all strings.
do i = 1, ras%num_strings
! This is equiavlent to F(J_{beta}) in Eq. 23 (except that we don't need to store
! the value for all strings at once).
alpha_beta_fac = 0.0_dp
! Get info for string_i.
sym_i = ras_strings(-1, i)
class_i = ras_strings(0, i)
ind_i = i - ras%cum_classes(class_i) - classes(class_i)%cum_sym(sym_i)
! Loop over all single excitations from string_i (-> string_j).
do excit_j = 1, ras_excit(i)%nexcit
full_ind_j = ras_excit(i)%excit_ind(excit_j)
class_j = ras_strings(0, full_ind_j)
par_1 = ras_excit(i)%par(excit_j)
ex2 = ras_excit(i)%orbs(:, excit_j)
BRR_ex2(1) = BRR(2 * ex2(1)) / 2
BRR_ex2(2) = BRR(2 * ex2(2)) / 2
alpha_beta_fac(full_ind_j) = alpha_beta_fac(full_ind_j) + &
par_1 * get_umat_el(BRR_ex2(2), BRR_ex1(2), BRR_ex2(1), BRR_ex1(1))
end do
! Loop over all classes allowed with class_i.
do j = 1, classes(class_i)%num_comb
class_j = classes(class_i)%allowed_combns(j)
sym_j = ieor(HFSym_ras, sym_i)
do ind_j = 1, classes(class_j)%num_sym(sym_j)
full_ind_j = ind_j + ras%cum_classes(class_j) + classes(class_j)%cum_sym(sym_j)
! If this is true then (kl) isn't a valid excitation from string j.
if (r(full_ind_j) == 0) cycle
v = 0.0_dp
! excited_class and excited_sym give the class and symmetry of the string
! which we get by exciting string j with (kl).
excited_class = excit_info(1, full_ind_j)
excited_sym = excit_info(2, full_ind_j)
! Loop over all classes connected to excited_class.
do m = 1, classes(excited_class)%num_comb
! Get the class and required symmetry.
class_m = classes(excited_class)%allowed_combns(m)
sym_m = ieor(HFSym_ras, excited_sym)
if (classes(class_m)%num_sym(sym_m) /= 0) then
min_ind = ras%cum_classes(class_m) + classes(class_m)%cum_sym(sym_m) + 1
max_ind = min_ind + classes(class_m)%num_sym(sym_m) - 1
v = v + dot_product(alpha_beta_fac(min_ind:max_ind), c(min_ind:max_ind, full_ind_j))
end if
end do
vec_out(class_j, class_i, sym_j)%elements(ind_j, ind_i) = &
vec_out(class_j, class_i, sym_j)%elements(ind_j, ind_i) + v
end do ! Over all strings in class_j with symmetry sym_j.
end do ! Over all classes allowed with class_i.
end do ! Over all strings.
end do ! Over all orbitals, l.
end do ! Over all orbitals, k.
! Deallocate arrays.
do class_i = 1, ras%num_classes
do sym_i = 0, 7
deallocate(factors(class_i, sym_i)%elements)
end do
end do
deallocate(alpha_beta_fac)
deallocate(c)
deallocate(r)
end subroutine perform_multiplication