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