subroutine generate_double_pcpp(nI, elec_map, ilut, nJ, excitMat, tParity, pGen)
implicit none
! given the initial determinant (both as nI and ilut), create a random double excitation
! given by nJ/ilutnJ/excitMat with probability pGen. tParity indicates the fermi sign
! picked up by applying the excitation operator
! Input: nI - determinant to excite from
! elec_map - map to translate electron picks to orbitals
! ilut - determinant to excite from in ilut format
! nJ - on return, excited determinant
! excitMat - on return, excitation matrix nI -> nJ
! tParity - on return, the parity of the excitation nI -> nJ
! pGen - on return, the probability of generating the excitation nI -> nJ
integer, intent(in) :: nI(nel)
integer, intent(in) :: elec_map(nel)
integer(n_int), intent(in) :: ilut(0:NIfTot)
integer, intent(out) :: nJ(nel)
integer, intent(out) :: ExcitMat(2, 2)
logical, intent(out) :: tParity
real(dp), intent(out) :: pGen
! temporary storage for the probabilities in each step
real(dp) :: pSGen1, pSGen2, pTGen1, pTGen2
! temporary storage for the probabilities in the swapped steps
real(dp) :: pSSwap1, pSSwap2, pTSwap1, pTSwap2
! temporary storage for the unmapped electrons
integer :: umElec1, umElec2, swapOrb1
! chosen source/target orbitals
integer :: src1, src2
integer :: tgt1, tgt2
integer :: elec1, elec2
! symmetry + spin to enforce for the last orbital
type(Symmetry) :: tgtSym
integer :: tgt1MS, tgt2MS
call double_elec_one_sampler%sample(umElec1, pSGen1)
src1 = elec_map(umElec1)
! in very rare cases, no mapping is possible in the first place
! then, abort
if (invalid_mapping(src1)) then
! the pgen here is just for bookkeeping purpose, it is never
! used
pGen = pSGen1
return
end if
call double_elec_two_sampler(src1)%sample(umElec2, pSGen2)
src2 = elec_map(umElec2)
! it is possible to not be able to map the second electron if
! the first mapping occupied the only available slot
if (invalid_mapping(src2, src1)) then
pGen = pSGen1 * pSGen2
return
end if
if (src2 < src1) call swap(src1, src2)
! we could also have drawn the electrons the other way around
pSSwap1 = double_elec_one_sampler%get_prob(umElec2)
swapOrb1 = elec_map(umElec2)
pSSwap2 = double_elec_two_sampler(swapOrb1)%get_prob(umElec1)
! we now have industuingishible src1/src2, add the probabilites
! for drawing them either way
pGen = pSGen1 * pSGen2 + pSSwap1 * pSSwap2
! decide which spins the tgt orbs shall have
! default to: src1 has the same spin as tgt1
tgt1MS = getSpinIndex(src1)
tgt2MS = getSpinIndex(src2)
! if we have opposite spins, chose their distribution randomly
if (G1(src1)%MS /= G1(src2)%MS) then
if (genrand_real2_dSFMT() < 0.5) call swap(tgt1MS, tgt2MS)
pGen = pGen * 0.5
end if
call double_hole_one_sampler(src1, tgt1MS)%sample(tgt1, pTGen1)
! update generation probability so far to ensure it has a valid value on return in any case
if (abort_excit(tgt1)) then
pGen = pGen * pTGen1
return
end if
! we need a specific symmetry now
tgtSym = getTgtSym(tgt1)
call double_hole_two_sampler(src2, tgtSym%s, tgt2MS)%sample(tgt2, pTGen2)
if (abort_excit(tgt2, tgt1)) then
pGen = pGen * pTGen1 * pTGen2
return
end if
! Update the generation probability
! We could have drawn the target orbitals the other way around
! -> adapt pGen
pTSwap1 = double_hole_one_sampler(src1, tgt2MS)%get_prob(tgt2)
tgtSym = getTgtSym(tgt2)
pTSwap2 = double_hole_two_sampler(src2, tgtSym%s, tgt1MS)%get_prob(tgt1)
pGen = pGen * (pTGen1 * pTGen2 + pTSwap1 * pTSwap2)
! generate the output determinant
elec1 = binary_search_first_ge(nI, src1)
elec2 = binary_search_first_ge(nI, src2)
call make_double(nI, nJ, elec1, elec2, tgt1, tgt2, ExcitMat, tParity)
contains
function getTgtSym(tgt) result(sym)
! return the symmetry of the last target orbital given a first
! target orbital tgt
! Input: tgt - first target orbital
! Output: sym - symmetry of the missing orbital
implicit none
integer, intent(in) :: tgt
type(symmetry) :: sym
sym = symprod(G1(src1)%Sym, G1(src2)%Sym)
sym = symprod(sym, symconj(G1(tgt)%Sym))
end function getTgtSym
function getTgtSpin(tgt) result(ms)
! return the spin of the last target orbital given a first target orbital
! Input: tgt - first target orbital
! Output: ms - spin index of the missing orbital:
! 0 - alpha
! 1 - beta
implicit none
integer, intent(in) :: tgt
integer :: ms
! if the electrons have the same spin, return the spin index of tgt
if (G1(src1)%MS == G1(src2)%MS) then
ms = getSpinIndex(tgt)
else
! else, the opposite spin index
ms = 1 - getSpinIndex(tgt)
end if
end function getTgtSpin
function invalid_mapping(src, src2) result(abort)
! check if the mapping was successful
! Input: src - electron we want to know about: did the mapping succeed?
! src2 - other chosen electron
! Output: abort - true if the mapping failed
implicit none
integer, intent(in) :: src
integer, optional, intent(in) :: src2
logical :: abort
abort = src == 0
if (present(src2)) abort = abort .or. (src == src2)
if (abort) then
nJ = 0
tParity = .false.
ExcitMat = 0
ExcitMat(1, 1) = src
if (present(src2)) ExcitMat(1, 2) = src2
end if
end function invalid_mapping
function abort_excit(tgt, tgt2) result(abort)
! check if the target orbital is valid
! Input: tgt - orbital we want to know about: Is an excitation to this possible
! tgt2 - second target orbital if already obtained
! Output: abort - true if there is no allowed excitation
implicit none
integer, intent(in) :: tgt
integer, optional, intent(in) :: tgt2
logical :: abort
abort = IsOcc(ilut, tgt) .or. (tgt == 0)
if (present(tgt2)) abort = abort .or. tgt == tgt2
if (abort) then
nJ = 0
ExcitMat(1, 1) = src1
ExcitMat(1, 2) = src2
ExcitMat(2, 1) = tgt
if (present(tgt2)) then
ExcitMat(2, 2) = tgt2
else
ExcitMat(2, 2) = 0
end if
tParity = .false.
end if
end function abort_excit
end subroutine generate_double_pcpp