subroutine InitFCIMC_CSF()
implicit none
integer(n_int), allocatable :: initSpace(:, :)
integer :: count, nUp, nOpen
integer :: i, j, lwork, proc
integer :: DetHash, pos, TotWalkersTmp, nI(nel), nJ(nel)
integer(n_int) :: ilutJ(0:NIfTot)
integer(n_int), allocatable :: openSubspace(:)
real(dp), allocatable :: S2(:, :), eigsImag(:), eigs(:), evs(:, :), void(:, :), work(:)
real(dp) :: normalization, rawWeight, HDiag, tmpSgn(lenof_sign)
HElement_t(dp) :: HOffDiag
integer :: err
real(dp) :: HFWeight(inum_runs)
logical :: tSuccess
character(*), parameter :: t_r = "InitFCIMC_CSF"
! get the number of open orbitals
nOpen = sum(nOccOrbs) - sum(nClosedOrbs)
! in a closed shell system, nothing to do
if (nOpen == 0) then
call InitFCIMC_HF()
return
end if
! first, set up the space considered for the CSF
call generateInitSpace()
if (allocated(openSubspace)) deallocate(openSubspace)
! we now have initSpace(:,:) with iluts belonging to all possible initial
! dets (i.e. all dets contributing to the target CSF) -> construct S2
allocate(S2(count, count))
do i = 1, count
do j = 1, count
S2(i, j) = S2Matel(initSpace(:, i), initSpace(:, j))
end do
end do
! prepare the diagonalization
allocate(eigs(count))
allocate(evs(count, count))
allocate(eigsImag(count))
allocate(work(1))
allocate(void(0, 0))
! workspace query, get how much tmp memory we need
call dgeev('N', 'V', count, S2, count, eigs, eigsImag, void, count, evs, count, work, -1, err)
! allocate work array
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
! diagonalize S2
call dgeev('N', 'V', count, S2, count, eigs, eigsImag, void, count, evs, count, work, lwork, err)
deallocate(void)
deallocate(work)
deallocate(eigsImag)
! transfer the eigenvector
do i = 1, count
if (abs(S2Init * (S2Init + 1) - eigs(i)) < eps) exit
end do
if (i > count) then
call stop_all(t_r, "Requested S2 eigenvalue does not exist")
end if
eigs = evs(:, i)
! normalization = minval(eigs)
rawWeight = sum(abs(eigs))
normalization = InitialPart / rawWeight
! refLoc = maxloc(eigs)
eigs = eigs * normalization
TotWalkers = 0
iStartFreeSlot = 1
iEndFreeSlot = 0
do i = 1, count
call decode_bit_det(nI, initSpace(:, i))
proc = DetermineDetNode(nel, nI, 0)
if (iProcIndex == proc) then
HDiag = get_diagonal_matel(nI, initSpace(:, i))
HOffDiag = get_off_diagonal_matel(nI, initSpace(:, i))
DetHash = FindWalkerHash(nI, size(HashIndex))
TotWalkersTmp = int(TotWalkers)
tmpSgn = eigs(i)
call encode_sign(initSpace(:, i), tmpSgn)
if (tHPHF) then
call FindDetSpinSym(nI, nJ, nel)
call encodebitdet(nJ, ilutJ)
! if initSpace(:,i) is not the det of the HPHF pair we are storing,
! skip this - the correct contribution will be stored once
! the spin-flipped version is stored
if (DetBitLT(initSpace(:, i), ilutJ, NIfD) == 1) cycle
end if
call AddNewHashDet(TotWalkersTmp, initSpace(:, i), DetHash, nI, HDiag, HOffDiag, pos, err)
TotWalkers = TotWalkersTmp
end if
! reset the reference?
end do
call hash_table_lookup(HFDet, ilutHF, nifd, HashIndex, CurrentDets, i, DetHash, tSuccess)
if (tSuccess) then
call extract_sign(CurrentDets(:, i), tmpSgn)
do i = 1, inum_runs
HFWeight(i) = mag_of_run(tmpSgn, i)
end do
else
HFWeight = 0.0_dp
end if
AllTotParts = InitialPart
AllTotPartsOld = InitialPart
OldAllAvWalkersCyc = InitialPart
OldAllHFCyc = HFWeight
OldAllNoatHF = HFWeight
! cleanup
deallocate(evs)
deallocate(eigs)
deallocate(S2)
deallocate(initSpace)
contains
!------------------------------------------------------------------------------------------!
subroutine generateOpenOrbIluts()
use IntegralsData, only: nfrozen
implicit none
count = 0
nUp = (nel + lms) / 2 - sum(nClosedOrbs) + nfrozen / 2
do i = 1, 2**nOpen - 1
if (popcnt(i) == nUp) then
count = count + 1
if (allocated(openSubspace)) openSubspace(count) = i
end if
end do
end subroutine generateOpenOrbIluts
!------------------------------------------------------------------------------------------!
subroutine generateInitSpace()
implicit none
integer :: nI(nel), nIBase(nel)
integer :: iEl, iElBase, iOpen
integer :: openOrbList(nel)
logical :: previousCont, nextCont, tClosed
character(*), parameter :: t_r = "generateInitSpace"
! create a list of all open-shell determinants with the correct spin+orbs
nUp = (nel + lms) / 2
call generateOpenOrbIluts()
allocate(openSubspace(count))
call generateOpenOrbIluts()
! convert open-shell-only iluts into full iluts
! use the reference to determine which orbitals shall participate
iElBase = 1
nIBase = 0
! generate a list of open orbitals
iOpen = 0
openOrbList = 0
do i = 1, nel
if (i == 1) then
previousCont = .false.
else
previousCont = FDet(i - 1) == FDet(i) - 1
end if
if (i == nel) then
nextCont = .false.
else
nextCont = FDet(i + 1) == FDet(i) + 1
end if
tClosed = .true.
! identify open orbitals, using the ordering: does the next one belong
! to the same spatial orb?
if (is_beta(FDet(i)) .and. .not. nextCont) then
iOpen = iOpen + 1
openOrbList(iOpen) = FDet(i)
tClosed = .false.
end if
! or the previous one?
if (is_alpha(FDet(i)) .and. .not. previousCont) then
iOpen = iOpen + 1
openOrbList(iOpen) = FDet(i)
tClosed = .false.
end if
! if the orbital is not open, it is closed
if (tClosed) then
nIBase(iElBase) = FDet(i)
iElBase = iElBase + 1
end if
end do
if (iOpen /= nOpen) then
write(stderr, *) "nOpen/iOpen conflict", FDet, openOrbList, iOpen, nOpen
call stop_all(t_r, "Error in determining open shell orbitals")
end if
allocate(initSpace(0:NIfTot, count))
! now, add the open-shell contribution
do i = 1, count
! start from the closed-shell base
nI = nIBase
iEl = iElBase
! for each open orb, add the electron
do j = 1, nOpen
! based on openSubspace(i), we are considering alpha/beta for this
! orb
if (btest(openSubspace(i), j - 1)) then
nI(iEl) = get_alpha(openOrbList(j))
else
nI(iEl) = get_beta(openOrbList(j))
end if
iEl = iEl + 1
end do
! encode the determinant to the initial space
call sort(nI)
call EncodeBitDet(nI, initSpace(:, i))
end do
end subroutine generateInitSpace
function S2Matel(ilutA, ilutB) result(matel)
implicit none
integer(n_int), intent(in) :: ilutA(0:NIfTot), ilutB(0:NIfTot)
integer(n_int) :: splus(0:NIfTot), sminus(0:NIfTot)
real(dp) :: matel
integer :: k, m, nI(nel), upOrb, downOrb
matel = 0.0_dp
if (DetBitEq(ilutA, ilutB, NIfD)) then
matel = matel + real(lms * (lms + 2), dp) / 4.0_dp
end if
! get the offdiag part of S2: S-S+
call decode_bit_det(nI, ilutA)
do k = 1, nel
! first, apply S- to all electrons (additively)
if (is_beta(nI(k))) then
sminus = ilutA
downOrb = get_alpha(nI(k))
clr_orb(sminus, nI(k))
set_orb(sminus, downOrb)
! now, apply S+ to all electrons (again, sum)
do m = 1, nel
! check if it yields 0
if (is_alpha(nI(m)) .or. m == k) then
splus = sminus
! if not, go on
if (m == k) then
clr_orb(splus, downOrb)
else
clr_orb(splus, nI(m))
end if
upOrb = get_beta(nI(m))
set_orb(splus, upOrb)
if (DetBitEq(splus, ilutB, NIFD)) matel = matel + 1.0_dp
end if
end do
end if
end do
end function S2Matel
end subroutine InitFCIMC_CSF