subroutine assign_reference_dets()
! Depending on the configuration we may have one, or multiple,
! reference determinants.
integer :: det(nel), orbs(nel), orb, orb2, norb
integer :: run, cc_idx, label_idx, i, j, found_orbs(inum_runs)
real(dp) :: energies(nel), hdiag
! for now guga only works with non-complex code
integer(n_int), allocatable :: excitations(:, :)
integer :: n_excits, ierr
real(dp), allocatable :: diag_energies(:)
logical, allocatable :: found_mask(:)
character(*), parameter :: this_routine = "assign_reference_dets"
! If the user has specified all of the (multiple) reference states,
! then just copy these across to the ilutRef array:
if (tMultipleInitialStates) then
tReplicaReferencesDiffer = .true.
do run = 1, inum_runs
ProjEDet(:, run) = initial_states(:, run)
call EncodeBitDet(ProjEDet(:, run), ilutRef(:, run))
end do
else if (tOrthogonaliseReplicas) then
tReplicaReferencesDiffer = .true.
! The first replica is just a normal FCIQMC simulation.
ilutRef(:, 1) = ilutHF
ProjEDet(:, 1) = HFDet
found_orbs = 0
! in the GUGA approach we have to change that simple
! single excitation of course or otherwise we get non-
! allowed CSF or the wrong STOT symmetry sector.
if (tGUGA) then
! run the exact single excitations on the HF det
! and find inum_runs lowest energetically ones..
! create all excitations from the HF
! do i remain in the same symmetry sector with this
! info??
! i think my exact hamil application routine does not
! deal with symmetry at all...
! i guess i have to change that for the real-space
! hubbard model implementation!
if (tHUB .or. tUEG .or. .not. (tNoBrillouin)) then
call actHamiltonian(ilutHF, CSF_Info_t(ilutHF), excitations, n_excits)
else
call actHamiltonian(ilutHF, CSF_Info_t(ilutHF), excitations, n_excits, .true.)
end if
! if no excitations possible... there is something wrong
if (.not. (n_excits > 0) .or. n_excits < inum_runs - 1) then
print *, "n_excits: ", n_excits
print *, "requested excited states:", inum_runs - 1
call write_guga_list(stdout, excitations(:, 1:n_excits))
call stop_all(this_routine, "not enough excitations from HF!")
end if
! and the choose the inum_runs - 1 lowest energetically
! excitations
! also have to calculate the correct diagonal element to
! compare energies
allocate(diag_energies(n_excits), stat=ierr)
allocate(found_mask(n_excits), stat=ierr)
found_mask = .true.
do i = 1, n_excits
diag_energies(i) = real(calcDiagMatEleGUGA_ilut(excitations(:, i)), dp)
end do
! can i sort the excitation list, according to energies?
do run = 2, inum_runs
! find the minimum energy
i = minloc(diag_energies, 1, found_mask)
! dont find that specific one again
found_mask(i) = .false.
! assign the reference determinant(transform between
! guga and actual iluts!)
call convert_ilut_toNECI(excitations(:, i), ilutRef(:, run))
! and calc. the nI representation
call decode_bit_det(ProjEDet(:, run), ilutRef(:, run))
! that should be it.. should the diagonal element also
! be already stored within the ilut? is that done here?
end do
else
do run = 2, inum_runs
! Now we want to find the lowest energy single excitation with
! the same symmetry as the reference site.
do i = 1, nel
! Find the excitations, and their energy
orb = HFDet(i)
cc_idx = ClassCountInd(orb)
label_idx = SymLabelCounts2(1, cc_idx)
norb = OrbClassCount(cc_idx)
! nb. sltcnd_0 does not depend on the ordering of the det,
! so we don't need to do any sorting here.
energies(i) = 9999999.9_dp
do j = 1, norb
orb2 = SymLabelList2(label_idx + j - 1)
if (.not. (any(orb2 == HFDet) .or. any(orb2 == found_orbs))) then
det = HFDet
det(i) = orb2
hdiag = real(sltcnd_excit(det, Excite_0_t()), dp)
if (hdiag < energies(i)) then
energies(i) = hdiag
orbs(i) = orb2
end if
end if
end do
end do
! Which of the electrons that is excited gives the lowest energy?
i = minloc(energies, 1)
found_orbs(run) = orbs(i)
! Construct that determinant, and set it as the reference.
ProjEDet(:, run) = HFDet
ProjEDet(i, run) = orbs(i)
call sort(ProjEDet(:, run))
call EncodeBitDet(ProjEDet(:, run), ilutRef(:, run))
end do
end if
else if (tPreCond) then
do run = 1, inum_runs
ilutRef(:, run) = ilutHF
ProjEDet(:, run) = HFDet
end do
! And make sure that the rest of the code knows this
tReplicaReferencesDiffer = .true.
else
! This is the normal case. All simultions are essentially doing
! the same thing...
do run = 1, inum_runs
ilutRef(:, run) = ilutHF
ProjEDet(:, run) = HFDet
end do
! And make sure that the rest of the code knows this
tReplicaReferencesDiffer = .false.
end if
write(stdout, *) 'Generated reference determinants:'
do run = 1, inum_runs
call write_det(stdout, ProjEDet(:, run), .false.)
hdiag = real(get_helement(ProjEDet(:, run), ProjEDet(:, run), 0), dp)
write(stdout, '(" E = ", f16.9)') hdiag
end do
end subroutine assign_reference_dets