subroutine InitFCIMC_MP1()
real(dp) :: TotMP1Weight, amp, MP2Energy, PartFac, rat, r, energy_contrib
HElement_t(dp) :: HDiagtemp, HOffDiagtemp
integer :: iExcits, exflag, Ex(2, maxExcit), nJ(NEl), DetIndex, iNode
integer :: iInit, DetHash, ExcitLevel, run, part_type
integer(n_int) :: iLutnJ(0:NIfTot)
real(dp) :: NoWalkers, temp_sign(lenof_sign)
logical :: tAllExcitsFound, tParity
type(ll_node), pointer :: TempNode
character(len=*), parameter :: this_routine = "InitFCIMC_MP1"
integer(n_int) :: ilutG(0:nifguga)
integer(n_int), allocatable :: excitations(:, :)
integer :: i
#ifdef CMPLX_
call stop_all(this_routine, "StartMP1 currently does not work with complex walkers")
#endif
if (tReadPops) call stop_all(this_routine, "StartMP1 cannot work with with ReadPops")
if (tStartSinglePart) call stop_all(this_routine, "StartMP1 cannot work with StartSinglePart")
if (tRestartHighPop) call stop_all(this_routine, "StartMP1 cannot with with dynamically restarting calculations")
write(stdout, *) "Initialising walkers proportional to the MP1 amplitudes..."
if (tHPHF) then
if (.not. TestClosedShellDet(iLutHF)) then
call stop_all(this_routine, "Cannot use HPHF with StartMP1 if your reference is open-shell")
end if
end if
if (tUEG) then
!Parallel N^2M implementation of MP2 for UEG
call CalcUEGMP2()
end if
!First, calculate the total weight - TotMP1Weight
mp2energy = 0.0_dp
TotMP1Weight = 1.0_dp
iExcits = 0
tAllExcitsFound = .false.
if (tUEG) then
exflag = 2
else
exflag = 3
end if
Ex(:, :) = 0
! for GUGA, this whole excitation generation is different!
if (tGUGA) then
! should finally do some general routine, which does all this
! below...
call convert_ilut_toGUGA(ilutHF, ilutG)
call actHamiltonian(ilutG, CSF_Info_t(ilutG), excitations, iExcits)
do i = 1, iExcits
call convert_ilut_toNECI(excitations(:, i), ilutnJ)
call decode_bit_det(nJ, iLutnJ)
call return_mp1_amp_and_mp2_energy(nJ, iLutnJ, Ex, tParity, amp, energy_contrib)
TotMP1Weight = TotMP1Weight + abs(amp)
MP2Energy = MP2Energy + energy_contrib
end do
else
do while (.true.)
call GenExcitations3(HFDet, iLutHF, nJ, exflag, Ex, tParity, tAllExcitsFound, .false.)
if (tAllExcitsFound) exit !All excits found
call EncodeBitDet(nJ, iLutnJ)
if (tHPHF) then
!Working in HPHF Space. Check whether determinant generated is an 'HPHF'
if (.not. IsAllowedHPHF(iLutnJ)) cycle
end if
iExcits = iExcits + 1
call return_mp1_amp_and_mp2_energy(nJ, iLutnJ, Ex, tParity, amp, energy_contrib)
TotMP1Weight = TotMP1Weight + abs(amp)
MP2Energy = MP2Energy + energy_contrib
end do
end if
if ((.not. tHPHF .and. .not. tGUGA) .and. (iExcits /= (nDoubles + nSingles))) then
write(stderr, *) nDoubles, nSingles, iExcits
call stop_all(this_routine, "Not all excitations accounted for in StartMP1")
end if
write(stdout, "(A,2G25.15)") "MP2 energy calculated: ", MP2Energy, MP2Energy + Hii
if ((InitialPart.isclose.1._dp) .or. (InitialPart >= (InitWalkers * nNodes) - 50)) then
!Here, all the walkers will be assigned to the MP1 wavefunction.
!InitialPart = 1 by default
write(stdout, "(A)") "All walkers specified in input will be distributed according to the MP1 wavefunction."
write(stdout, "(A)") "Shift will be allowed to vary from the beginning"
write(stdout, "(A)") "Setting initial shift to equal MP2 correlation energy"
DiagSft = MP2Energy
!PartFac is the number of walkers that should reside on the HF determinant
!in an intermediate normalised MP1 wavefunction.
PartFac = (real(InitWalkers, dp) * real(nNodes, dp)) / TotMP1Weight
else
!Here, not all walkers allowed will be initialised to the MP1 wavefunction.
write(stdout, "(A,G15.5,A)") "Initialising ", InitialPart, " walkers according to the MP1 distribution."
write(stdout, "(A,G15.5)") "Shift will remain fixed until the walker population reaches ", InitWalkers * nNodes
!PartFac is the number of walkers that should reside on the HF determinant
!in an intermediate normalised MP1 wavefunction.
PartFac = real(InitialPart, dp) / TotMP1Weight
tSinglePartPhase(:) = .true.
end if
!Now generate all excitations again, creating the required number of walkers on each one.
! puh... for GUGA this gets messy to change below.. damn
DetIndex = 1
TotParts = 0.0
tAllExcitsFound = .false.
if (tUEG) then
exflag = 2
else
exflag = 3
end if
Ex(:, :) = 0
! figure out if the HF det gets store in the excitation list too?
! if yes I have to modify that all a bit, and maybe also in
! other parts of the NECI code ... todo
if (tGUGA) call stop_all(this_routine, "deprecated option with GUGA!")
do while (.true.)
call GenExcitations3(HFDet, iLutHF, nJ, exflag, Ex, tParity, tAllExcitsFound, .false.)
if (tAllExcitsFound) exit !All excits found
call EncodeBitDet(nJ, iLutnJ)
if (tHPHF) then
!Working in HPHF Space. Check whether determinant generated is an 'HPHF'
if (.not. IsAllowedHPHF(iLutnJ)) cycle
end if
iNode = DetermineDetNode(nel, nJ, 0)
if (iProcIndex == iNode) then
call return_mp1_amp_and_mp2_energy(nJ, iLutnJ, Ex, tParity, amp, energy_contrib)
amp = amp * PartFac
if (tRealCoeffByExcitLevel) ExcitLevel = FindBitExcitLevel(iLutnJ, iLutRef(:, 1), nEl)
if (tAllRealCoeff .or. &
& (tRealCoeffByExcitLevel .and. (ExcitLevel <= RealCoeffExcitThresh))) then
NoWalkers = amp
else
NoWalkers = int(amp)
rat = amp - real(NoWalkers, dp)
r = genrand_real2_dSFMT()
if (abs(rat) > r) then
if (amp < 0.0_dp) then
NoWalkers = NoWalkers - 1
else
NoWalkers = NoWalkers + 1
end if
end if
end if
if (abs(NoWalkers) > 1.0e-12_dp) then
call encode_det(CurrentDets(:, DetIndex), iLutnJ)
call clear_all_flags(CurrentDets(:, DetIndex))
do run = 1, inum_runs
temp_sign(run) = NoWalkers
end do
call encode_sign(CurrentDets(:, DetIndex), temp_sign)
! Store the diagonal matrix elements
HDiagtemp = get_diagonal_matel(nJ, iLutnJ)
HOffDiagtemp = get_off_diagonal_matel(nJ, iLutnJ)
call set_det_diagH(DetIndex, real(HDiagtemp, dp) - Hii)
call set_det_offdiagH(DetIndex, HOffDiagtemp)
if (associated(lookup_supergroup_indexer)) then
call set_supergroup_idx(DetIndex, lookup_supergroup_indexer%idx_nI(nJ))
end if
! store the determinant
call store_decoding(DetIndex, nJ)
if (tTruncInitiator) then
!Set initiator flag if needed (always for HF)
call CalcParentFlag(DetIndex, iInit)
end if
DetHash = FindWalkerHash(nJ, nWalkerHashes)
TempNode => HashIndex(DetHash)
! If the first element in the list has not been used.
if (TempNode%Ind == 0) then
TempNode%Ind = DetIndex
else
do while (associated(TempNode%Next))
TempNode => TempNode%Next
end do
allocate(TempNode%Next)
nullify (TempNode%Next%Next)
TempNode%Next%Ind = DetIndex
end if
nullify (TempNode)
DetIndex = DetIndex + 1
do part_type = 1, lenof_sign
TotParts(part_type) = TotParts(part_type) + abs(NoWalkers)
end do
end if
end if !End if desired node
end do
!Now for the walkers on the HF det
if (iRefProc(1) == iProcIndex) then
! dont have to change this below, should also work for the GUGAc
if (tAllRealCoeff .or. tRealCoeffByExcitLevel) then
NoWalkers = PartFac
else
NoWalkers = int(PartFac)
rat = PartFac - real(NoWalkers, dp)
if (rat < 0.0_dp) &
call stop_all(this_routine, "Should not have negative weight on HF")
r = genrand_real2_dSFMT()
if (abs(rat) > r) NoWalkers = NoWalkers + 1
end if
if (abs(NoWalkers) > 1.0e-12_dp) then
call encode_det(CurrentDets(:, DetIndex), iLutHF)
call clear_all_flags(CurrentDets(:, DetIndex))
do run = 1, inum_runs
temp_sign(run) = NoWalkers
end do
call encode_sign(CurrentDets(:, DetIndex), temp_sign)
if (tTruncInitiator) then
!Set initiator flag (always for HF)
do run = 1, inum_runs
call set_flag(CurrentDets(:, DetIndex), get_initiator_flag(run))
end do
end if
call set_det_diagH(DetIndex, 0.0_dp)
call set_det_offdiagH(DetIndex, h_cast(0_dp))
! store the determinant
call store_decoding(DetIndex, HFDet)
! Now add the Hartree-Fock determinant (not with index 1).
DetHash = FindWalkerHash(HFDet, nWalkerHashes)
TempNode => HashIndex(DetHash)
! If the first element in the list has not been used.
if (TempNode%Ind == 0) then
TempNode%Ind = DetIndex
else
do while (associated(TempNode%Next))
TempNode => TempNode%Next
end do
allocate(TempNode%Next)
nullify (TempNode%Next%Next)
TempNode%Next%Ind = DetIndex
end if
nullify (TempNode)
DetIndex = DetIndex + 1
do run = 1, inum_runs
TotParts(run) = TotParts(run) + abs(NoWalkers)
NoatHF(run) = NoWalkers
end do
else
call stop_all(this_routine, "No walkers initialised on the HF det with StartMP1")
end if
else
NoatHF(:) = 0.0_dp
end if
TotWalkers = DetIndex - 1 !This is the number of occupied determinants on each node
TotWalkersOld = TotWalkers
!Set local&global variables
TotPartsOld = TotParts
call mpisumall(TotParts, AllTotParts)
call mpisumall(NoatHF, AllNoatHF)
call mpisumall(TotWalkers, AllTotWalkers)
OldAllNoatHF = AllNoatHF
do run = 1, inum_runs
OldAllHFCyc(run) = AllNoatHF(run)
OldAllAvWalkersCyc(run) = AllTotParts(run)
end do
AllTotWalkersOld = AllTotWalkers
AllTotPartsOld = AllTotParts
iter_data_fciqmc%tot_parts_old = AllTotPartsOld
AllNoAbortedOld = 0.0_dp
end subroutine InitFCIMC_MP1