SUBROUTINE CalcApproxpDoubles()
implicit none
real(dp) :: denom
INTEGER :: iTotal
integer :: nSingles, nDoubles, nSing_spindiff1, nDoub_spindiff1, nDoub_spindiff2
integer :: nTot
character(*), parameter :: this_routine = "CalcApproxpDoubles"
integer(n_int), allocatable :: dummy_list(:, :)
nSingles = 0
nDoubles = 0
if (tReltvy) then
nSing_spindiff1 = 0
nDoub_spindiff1 = 0
nDoub_spindiff2 = 0
pDoub_spindiff1 = 0.0_dp
pDoub_spindiff2 = 0.0_dp
end if
!NSing=Number singles from HF, nDoub=No Doubles from HF
write(stdout, "(A)") " Calculating approximate pDoubles for use with &
&excitation generator by looking a excitations from &
&reference."
exflag = 3
if (tReltvy) then
write(stdout, *) "Counting magnetic excitations"
! subroutine CountExcitations4(nI, minRank, maxRank, minSpinDiff, maxSpinDiff, tot)
call CountExcitations4(hfdet, 1, 1, 0, 0, nSingles)
call CountExcitations4(hfdet, 1, 1, 1, 1, nSing_spindiff1)
call CountExcitations4(hfdet, 2, 2, 0, 0, nDoubles)
call CountExcitations4(hfdet, 2, 2, 1, 1, nDoub_spindiff1)
call CountExcitations4(hfdet, 2, 2, 2, 2, nDoub_spindiff2)
call CountExcitations4(hfdet, 1, 2, 0, 2, nTot)
ASSERT(nTot == (nSingles + nSing_spindiff1 + nDoubles + nDoub_spindiff1 + nDoub_spindiff2))
iTotal = nSingles + nDoubles + nSing_spindiff1 + nDoub_spindiff1 + nDoub_spindiff2
else
if (tKPntSym) THEN
if (t_k_space_hubbard) then
! change this to the new implementation
call gen_all_excits_k_space_hubbard(HFDet, nDoubles, dummy_list)
else
call enumerate_sing_doub_kpnt(exFlag, .false., nSingles, nDoubles, .false.)
end if
else
call CountExcitations3(hfdet, exflag, nSingles, nDoubles)
end if
iTotal = nSingles + nDoubles
end if
IF (tHub .or. tUEG) THEN
IF (tReal) THEN
write(stdout, *) "Since we are using a real-space hubbard model, only single excitations are connected &
& and will be generated."
pDoubles = 0.0_dp
if (tReltvy) then
pDoub_spindiff1 = 0.0_dp
pDoub_spindiff2 = 0.0_dp
pSingles = real(nSingles, dp) / real(nSingles + nSing_spindiff1, dp)
pSing_spindiff1 = 1.0_dp - pSingles
else
pSingles = 1.0_dp
end if
return
ELSE
write(stdout, *) "Since we are using a momentum-space hubbard model/UEG, only double excitaitons &
& are connected and will be generated."
pSingles = 0.0_dp
if (tReltvy) then
pSing_spindiff1 = 0.0_dp
pDoubles = real(nDoubles, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
pDoub_spindiff1 = real(nDoub_spindiff1, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
pDoub_spindiff2 = 1.0_dp - pDoubles - pDoub_spindiff1
else
pDoubles = 1.0_dp
end if
return
end if
else if (tNoSingExcits) then
pSingles = 0.0_dp
if (tReltvy) then
pSing_spindiff1 = 0.0_dp
pDoubles = real(nDoubles, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
pDoub_spindiff1 = real(nDoub_spindiff1, dp) / real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
pDoub_spindiff2 = 1.0_dp - pDoubles - pDoub_spindiff1
else
pDoubles = 1.0_dp
end if
write(stdout, *) "Only double excitations will be generated"
return
end if
write(stdout, "(I7,A,I7,A)") nDoubles, " double excitations, and ", nSingles, &
" single excitations found from reference. This will be used to calculate pDoubles."
IF (abs(SinglesBias - 1.0_dp) > 1.0e-12_dp) THEN
write(stdout, *) "Singles Bias detected. Multiplying single excitation connectivity of HF determinant by ", &
SinglesBias, " to determine pDoubles."
end if
IF ((nSingles == 0) .or. (nDoubles == 0)) THEN
write(stdout, *) "Number of singles or doubles found equals zero. pDoubles will be set to 0.95. Is this correct?"
pDoubles = 0.95_dp
pSingles = 0.05_dp
return
else if ((nSingles < 0) .or. (nDoubles < 0)) then
call stop_all("CalcApproxpDoubles", &
"Number of singles, doubles or Yamanouchi symbols &
&found to be a negative number. Error here.")
end if
if (tReltvy) then
denom = real(nSingles + nSing_spindiff1, dp) * SinglesBias &
+ real(nDoubles + nDoub_spindiff1 + nDoub_spindiff2, dp)
pSingles = real(nSingles, dp) * SinglesBias / denom
pSing_spindiff1 = real(nSing_spindiff1, dp) * SinglesBias / denom
pDoubles = real(nDoubles, dp) / denom
pDoub_spindiff1 = real(nDoub_spindiff1, dp) / denom
pDoub_spindiff2 = 1.0_dp - pSingles - pSing_spindiff1 &
- pDoubles - pDoub_spindiff1 - pDoub_spindiff2
else
denom = real(nSingles, dp) * SinglesBias + real(nDoubles, dp)
pSingles = real(nSingles, dp) * SinglesBias / denom
pDoubles = 1.0_dp - pSingles
end if
IF (abs(SinglesBias - 1.0_dp) > 1.0e-12_dp) THEN
write (stdout, '("pDoubles set to ", f14.6, &
&" rather than (without bias): ", f14.6)') &
pDoubles, real(nDoubles, dp) / real(iTotal, dp)
write (stdout, '("pSingles set to ", f14.6, &
&" rather than (without bias): ", f14.6)') &
pSingles, real(nSingles, dp) / real(iTotal, dp)
! write(stdout,"(A,F14.6,A,F14.6)") "pDoubles set to: ",pDoubles, " rather than (without bias): ", &
! & real(nDoub,dp)/real(iTotal,dp)
ELSE
if (tReltvy) then
write(stdout, '(A)') " Where s and t are alpha or beta spin function labels: "
write(stdout, '(A30,F14.6)') " pSingles(s->s) set to: ", pSingles
write(stdout, '(A30,F14.6)') " pSingles(s->s') set to: ", pSing_spindiff1
write(stdout, '(A30,F14.6)') " pDoubles(st->st) set to: ", pDoubles
write(stdout, '(A30,F14.6)') " pDoubles(st->s't) set to: ", pDoub_spindiff1
write(stdout, '(A30,F14.6)') " pDoubles(st->s't') set to: ", pDoub_spindiff2
else
write(stdout, '(A,F14.6)') " pDoubles set to: ", pDoubles
write(stdout, '(A,F14.6)') " pSingles set to: ", pSingles
end if
end if
if (allocated(pSinglesIn) .or. allocated(pDoublesIn) .or. allocated(pTriplesIn)) then
if (allocated(pSinglesIn) .and. allocated(pDoublesIn)) then
call stop_all(this_routine, 'It is not possible to define pSingles and pDoubles.')
else if (.not. (allocated(pSinglesIn) .or. allocated(pDoublesIn))) then
call stop_all(this_routine, 'One of pSingles or pDoubles is required.')
end if
if (t_mol_3_body) then
! We allow the users to input absolute values for the probabilities.
! Note that pSingles and pDoubles are internally conditional probabilities.
! Even for triple we still have that `pSingles + pDoubles .isclose. 1.0`.
! It first decides upon triple excitation or something else and then about singles or doubles.
! We have to convert the absolute probabilities into conditional ones.
if (.not. allocated(pTriplesIn)) then
call stop_all(this_routine, "pTriples is required as input.")
else
pTriples = pTriplesIn
if (allocated(pSinglesIn)) then
if (pTriples + pSinglesIn > 1.0_dp) call stop_all(this_routine, "pTriplesIn + pSinglesIn > 1.0_dp")
pSingles = pSinglesIn / (1.0_dp - pTriplesIn)
pDoubles = 1.0_dp - pSingles
write(stdout, '(" Using the input value of pSingles:",1x, f14.6)') pSinglesIn
else if (allocated(pDoublesIn)) then
if (pTriples + pDoublesIn > 1.0_dp) call stop_all(this_routine, "pTriplesIn + pDoublesIn > 1.0_dp")
pDoubles = pDoublesIn / (1.0_dp - pTriplesIn)
pSingles = 1.0_dp - pDoubles
write(stdout, '(" Using the input value of pDoubles:",1x, f14.6)') pDoublesIn
end if
end if
else
if (allocated(pTriplesIn)) then
call stop_all(this_routine, "pTriples can only be given if triple excitations are performed.")
else if (allocated(pSinglesIn)) then
pSingles = pSinglesIn
pDoubles = 1.0_dp - pSingles
write(stdout, '(" Using the input value of pSingles:",1x, f14.6)') pSingles
else if (allocated(pDoublesIn)) then
pDoubles = pDoublesIn
pSingles = 1.0_dp - pDoubles
write(stdout, '(" Using the input value of pDoubles:",1x, f14.6)') pDoubles
end if
end if
end if
if (allocated(pParallelIn)) then
write(stdout, '(" Using the input value of pParallel:",1x, f14.6)') pParallelIn
pParallel = pParallelIn
end if
END SUBROUTINE CalcApproxpDoubles