subroutine create_ham_diag_direct_ci(ras, classes, ras_strings, ham_diag)
! The Davidson routine needs the Hamiltonian diagonal for the preconditioner.
! This routine creates and stores this array in the RAS space being used.
use Determinants, only: get_helement
use SystemData, only: nel
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)
real(dp), intent(inout) :: ham_diag(:)
integer :: class_i_ind, class_j_ind, sym_i_ind, sym_j_ind
integer :: class_i, class_j, j, sym_i, sym_j
integer :: ind_i, ind_j, full_ind_i, full_ind_j
integer :: counter, k
integer :: string_i(tot_nelec), string_j(tot_nelec), nI(nel)
counter = 0
! Over all classes.
do class_i = 1, ras%num_classes
class_i_ind = ras%cum_classes(class_i)
! Over all classes connected to class_i.
do j = 1, classes(class_i)%num_comb
class_j = classes(class_i)%allowed_combns(j)
class_j_ind = ras%cum_classes(class_j)
! Over all symmetry labels.
do sym_i = 0, 7
! The symmetry label for string_j is fixed by the label for string_i.
sym_j = ieor(HFSym_ras, sym_i)
! The index of the first string with this symmetry and class.
sym_i_ind = class_i_ind + classes(class_i)%cum_sym(sym_i)
sym_j_ind = class_j_ind + classes(class_j)%cum_sym(sym_j)
! If there are no strings with this symmetry and class.
if (classes(class_i)%num_sym(sym_i) == 0) cycle
if (classes(class_j)%num_sym(sym_j) == 0) cycle
! Over all strings with this symmetry and class, for string_i.
do ind_i = 1, classes(class_i)%num_sym(sym_i)
full_ind_i = sym_i_ind + ind_i
! Lookup the string corresponding to this address.
string_i = ras_strings(1:tot_nelec, full_ind_i)
! Over all strings with this symmetry and class, for string_j.
do ind_j = 1, classes(class_j)%num_sym(sym_j)
counter = counter + 1
full_ind_j = sym_j_ind + ind_j
string_j = ras_strings(1:tot_nelec, full_ind_j)
! Beta string.
nI(1:tot_nelec) = string_i * 2 - 1
! Alpha string.
nI(tot_nelec + 1:nel) = string_j * 2
! Replace all orbital numbers, orb, with the true orbital
! numbers, BRR(orb). Also, sort this list.
do k = 1, nel
nI(k) = BRR(nI(k))
end do
call sort(nI)
! Finally calculate and store the corresponding element.
ham_diag(counter) = get_helement(nI, nI, 0)
end do
end do
end do
end do
end do
end subroutine create_ham_diag_direct_ci