#include "macros.h" MODULE SymExcit3 ! This module contains excitation generators able to enumerate all possible excitations given a starting determinant. ! Unlike symexcit.F90 however, these excitation generators are able to deal with cases where the alpha and beta orbitals ! have different symmetries. This is particularly relevant when dealing with certain unrestricted cases, or when we ! are truncating (or freezing) orbitals in such a way as to remove different alpha symm irreps from the beta. use SystemData, only: NEl, G1, nBasis, tNoSymGenRandExcits use bit_rep_data, only: NIfTot use constants, only: n_int, maxExcit, stdout USE GenRandSymExcitNUMod, only: SymLabelList2, SymLabelCounts2, ClassCountInd, ScratchSize use SymExcitDataMod, only: SpinOrbSymLabel use get_excit, only: make_double use sort_mod, only: sort use util_mod, only: operator(.implies.), stop_all use DetBitOps, only: ilut_lt, ilut_gt, EncodeBitDet use excit_mod, only: FindExcitDet IMPLICIT NONE CONTAINS SUBROUTINE CountExcitations3(nI, exflag, nSingleExcits, nDoubleExcits) ! This routine simply counts the excitations in terms of single and doubles from the nI determinant. ! The exflag sent through indicates which should be counted - exflag=1 means only singles, exflag=2 means ! only doubles, and anything else both are counted. USE SymData, only: nSymLabels USE SystemData, only: ElecPairs, tFixLz, iMaxLz USE GenRandSymExcitNUMod, only: PickElecPair, construct_class_counts, ClassCountInd, ScratchSize INTEGER :: nSingleExcits, nDoubleExcits, Symi, i, Spini, nI(NEl) INTEGER :: iSpn, Elec1Ind, Elec2Ind, SymProduct, exflag INTEGER :: Syma, Symb, Spina, Spinb, StartSpin, EndSpin INTEGER :: ClassCount2(ScratchSize), SumMl INTEGER :: ClassCountUnocc2(ScratchSize) INTEGER :: StartMl, EndMl, Mla, Mlb CALL construct_class_counts(nI, ClassCount2, ClassCountUnocc2) ! This sets up arrays containing the number of occupied and unoccupied in each symmetry. ! ClassCounts2(1,:)=No alpha occupied, ClassCounts2(2,:)=No Beta occupied. ! ClassCountsUnocc2(1,:)=No alpha unocc, ClassCounts2Unocc2(2,:)=No Beta unocc. ! The second index of these arrays referrs to the symmetry (0 -> 7). ! Only counting. Run through each occupied orbital, and count the number of spin and symmetry allowed orbitals it ! may be excited to. nSingleExcits = 0 nDoubleExcits = 0 IF (exflag /= 2) THEN ! Count the singles. ! Take each electron and find out the number of symmetry allowed orbitals it may be excited to. do i = 1, NEl Symi = SpinOrbSymLabel(nI(i)) IF ((G1(nI(i))%Ms) == -1) Spini = 2 ! G1(i)%Ms is -1 for beta, and 1 for alpha. IF ((G1(nI(i))%Ms) == 1) Spini = 1 ! Translate this into 1 for alpha and 2 for beta ! for the ClassCount arrays. IF (tFixLz) THEN Mla = G1(nI(i))%Ml ELSE Mla = 0 end if ! This electron in orbital of SymI and SpinI can only be excited to orbitals with the same spin and symmetry. ! Then add in the number of unoccupied orbitals with the same spin and symmetry to which each electron may be excited. nSingleExcits = nSingleExcits + ClassCountUnocc2(ClassCountInd(Spini, Symi, Mla)) end do end if ! This is the end of the singles. ! write(stdout,*) 'Number of singles',nSingleExcits ! For the doubles, first pick an electron pair i,j. ! Based on these orbitals, run through each spin and each symmetry - take this to be orbital a. ! Multiply the number with these symmetries by the number of possible b orbitals which correspond. ! Do this for all a and then all i,j pairs. IF (exflag /= 1) THEN do i = 1, ElecPairs ! iSpn=2 for alpha beta pair, ispn=3 for alpha alpha pair and ispn=1 for beta beta pair. CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, i) StartSpin = 1 EndSpin = 2 IF (iSpn == 3) EndSpin = 1 IF (iSpn == 1) StartSpin = 2 do Spina = StartSpin, EndSpin ! Run through both spins, orbital a may be alpha or beta. IF (iSpn == 2) THEN ! Spin of orbital b should be opposite to orbital a. IF (Spina == 1) Spinb = 2 IF (Spina == 2) Spinb = 1 ELSE ! Spin of orbital b should be the same as orbital a. IF (Spina == 1) Spinb = 1 IF (Spina == 2) Spinb = 2 end if do Syma = 0, nSymLabels - 1 ! Need to work out the symmetry of b, given the symmetry of a (Sym). Symb = IEOR(Syma, SymProduct) IF (tFixLz) THEN StartMl = -iMaxLz EndMl = iMaxLz ELSE StartMl = 0 EndMl = 0 end if do Mla = StartMl, EndMl Mlb = SumMl - Mla !Will be 0 if no Lz, otherwise we need Mla + Mlb = Mli + Mlj = SumMl IF (ABS(Mlb) <= iMaxLz) THEN IF ((Spina == Spinb) .and. (Syma == Symb) .and. (Mla == Mlb)) THEN ! If the spin and spatial symmetries of a and b are the same ! there will exist a case where Orba = Orbb, want to remove this. nDoubleExcits = nDoubleExcits + (ClassCountUnocc2(ClassCountInd(Spina, Syma, Mla)) & * (ClassCountUnocc2(ClassCountInd(Spinb, Symb, Mlb)) - 1)) ELSE nDoubleExcits = nDoubleExcits + (ClassCountUnocc2(ClassCountInd(Spina, Syma, Mla)) & * ClassCountUnocc2(ClassCountInd(Spinb, Symb, Mlb))) end if end if end do end do end do end do nDoubleExcits = nDoubleExcits / 2 end if ENDSUBROUTINE CountExcitations3 SUBROUTINE GenExcitations3(nI, iLut, nJ, exflag, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only) ! This routine finds in turn, every possible excitation from determinant nI. ! The excited determinant is then returned as nJ. ! exflag indicates which excitations we want to find. exflag=1 - only singles are returned, exflag=2 - only ! doubles are returned and anything else returns the singles followed by the doubles. ! ExcitMat3 holds the orbitals involved in the excitation. ! If an excitation matrix of 0's is passed through, the first single or double is found. ! After this, the routine reads in the ExcitMat and finds the next excitation after this. ! ExcitMat(1,*) are the orbitals in the determinant to vacate from nI (the i,j pair) ! ExcitMat(2,*) are the orbitals to occupy in nJ (the a,b pair) (not the index, but the actual orbital) ! If tParity is true, two orbitals need to be switched in order to better represent the excitation, therefore a ! negative sign must be included when finding the H element. ! When there are no more symmetry allowed excitations, tAllExcitFound becomes true. INTEGER(KIND=n_int), intent(in) :: iLut(0:NIfTot) INTEGER, intent(in) :: nI(NEl) integer, intent(out) :: nJ(NEl) integer, intent(inout) :: ExcitMat3(2, 2), exflag LOGICAL, intent(out) :: tAllExcitFound, tParity LOGICAL, intent(in) :: ti_lt_a_only tAllExcitFound = .false. IF (exflag == 2) THEN ! Just generate doubles CALL GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only) ELSE ! Generate singles, returning Orbi and Orba as non-zero, but keeping the others 0. CALL GenSingleExcit(nI, iLut, nJ, exflag, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only) ! When the last single is input, providing exflag is not 1, the first double is then found ! and from then on GenDoubleExcit is called. end if ENDSUBROUTINE GenExcitations3 SUBROUTINE GenSingleExcit(nI, iLut, nJ, exflag, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only) ! Despite being fed four indices, this routine finds single excitations. Orbi -> Orba. (Orbj and Orbb remain 0). ! Feeding in 0 indices indicates it is the first excitation that needs to be found. ! The single excitation goes from orbital i to a, from determinant nI to nJ. ! When the last single is found it then finds the first double excitation, unless exflag=1 in which tAllExcitFound ! becomes true and no more excitations are generated. use SystemData, only: tFixLz use constants, only: bits_n_int INTEGER :: nI(NEl), Orbi, Orba, Symi, nJ(NEl) INTEGER(KIND=n_int) :: iLut(0:NIfTot) INTEGER :: NoOcc, ExcitMat3(2, 2), exflag, SymInd, Spina, Mla LOGICAL :: tInitOrbsFound, tParity, tAllExcitFound, tEndaOrbs, ti_lt_a_only, tAux INTEGER, SAVE :: OrbiIndex, OrbaIndex, Spini, NewSym, Mli tInitOrbsFound = .false. Orbi = ExcitMat3(1, 1) Orba = ExcitMat3(2, 1) IF ((Orbi == 0) .or. (Orba == 0)) THEN ! Want to find the first excitation. OrbiIndex = 1 Orbi = nI(OrbiIndex) ! Take the first occupied orbital Symi = SpinOrbSymLabel(Orbi) ! and find its spin and spat symmetries. IF ((G1(Orbi)%Ms) == -1) Spini = 2 IF ((G1(Orbi)%Ms) == 1) Spini = 1 IF (tFixLz) THEN Mli = G1(Orbi)%Ml ELSE Mli = 0 end if OrbaIndex = SymLabelCounts2(1, ClassCountInd(Spini, Symi, Mli)) ! Start considering a at the first allowed symmetry. ELSE Orbi = nI(OrbiIndex) ! Begin by using the same i as last time - check if there are any ! more possible excitations from this. ! At this stage, OrbaIndex is the a from the previous excitation. SymInd = ClassCountInd(Spini, SpinOrbSymLabel(Orbi), Mli) IF (OrbaIndex == (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) THEN !Orba was the last in the symmetry block. Do not allow OrbaIndex+1 ! Either we're got to the final spin symmetry, or the next orbital after !Orba does not have the same symmetry as Orbi. ! Need to move onto the next i, and find a new a to match. OrbiIndex = OrbiIndex + 1 IF (OrbiIndex <= NEl) THEN Orbi = nI(OrbiIndex) Symi = SpinOrbSymLabel(Orbi) IF ((G1(Orbi)%Ms) == -1) Spini = 2 IF ((G1(Orbi)%Ms) == 1) Spini = 1 IF (tFixLz) THEN Mli = G1(Orbi)%Ml ELSE Mli = 0 end if OrbaIndex = SymLabelCounts2(1, ClassCountInd(Spini, Symi, Mli)) ELSE IF (exflag /= 1) THEN ExcitMat3(:, :) = 0 CALL GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only) exflag = 2 ELSE tAllExcitFound = .true. tInitOrbsFound = .true. end if end if ELSE ! There are more possible excitations from orbital a, simply check the next orbital after the current a. OrbaIndex = OrbaIndex + 1 Symi = SpinOrbSymLabel(Orbi) end if end if do while (.not. tInitOrbsFound) tEndaOrbs = .false. IF (OrbiIndex > NEl) THEN ! If we've read in the last single, set orbi, orbj, orba, and orbb to 0 and call gendoubleexcit. IF (exflag /= 1) THEN ExcitMat3(:, :) = 0 CALL GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, ti_lt_a_only) exflag = 2 ELSE tAllExcitFound = .true. end if EXIT end if ! To find Orba, take the first in SymLabelList2 with the same symmetry and spin. ! SymLabelCounts2(spin,1,symmetry) gives the index in SymLabelList2 where that spin and symmetry starts. IF (OrbaIndex > nBasis) THEN tEndaOrbs = .true. ELSE tEndaOrbs = .false. Orba = SymLabelList2(OrbaIndex) end if SymInd = ClassCountInd(Spini, SpinOrbSymLabel(Orbi), Mli) ! Need to also make sure orbital a is unoccupied, so make sure the orbital is not in nI. NoOcc = 0 IF (.not. tEndaOrbs) THEN do while ((BTEST(iLut((Orba - 1) / bits_n_int), MOD((Orba - 1), bits_n_int))) .or. & (ti_lt_a_only .and. (Orba < Orbi))) ! While this is true, Orba is occupied, so keep incrementing Orba until it is not. NoOcc = NoOcc + 1 IF (OrbaIndex + NoOcc > nBasis) THEN !We have reached the end of all a orbitals. Now we need to pick a new i tEndaOrbs = .true. EXIT ELSE Orba = SymLabelList2(OrbaIndex + NoOcc) IF ((OrbaIndex + NoOcc) > (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) EXIT end if end do end if tAux = .false. IF (.not. tEndaOrbs) THEN ! Then check we have not overrun the symmetry block while skipping the occupied orbitals. NewSym = SpinOrbSymLabel(Orba) IF ((G1(Orba)%Ms) == -1) Spina = 2 IF ((G1(Orba)%Ms) == 1) Spina = 1 IF (tFixLz) THEN Mla = G1(Orba)%Ml ELSE Mla = 0 end if IF (NewSym == Symi .and. (Spina == Spini) .and. (Mli == Mla)) THEN ! If not, then these are the new Orbi and Orba. tInitOrbsFound = .true. OrbaIndex = OrbaIndex + NoOcc ELSE tAux = .true. end if end if ! If we have, move onto the next occupied orbital i, no symmetry allowed single excitations exist from the first. IF (tAux .or. tEndaOrbs) THEN OrbiIndex = OrbiIndex + 1 IF (OrbiIndex <= NEl) THEN Orbi = nI(OrbiIndex) Symi = SpinOrbSymLabel(Orbi) ! and find its spin and spat symmetries. IF ((G1(Orbi)%Ms) == -1) Spini = 2 IF ((G1(Orbi)%Ms) == 1) Spini = 1 IF (tFixLz) THEN Mli = G1(Orbi)%Ml ELSE Mli = 0 end if OrbaIndex = SymLabelCounts2(1, ClassCountInd(Spini, Symi, Mli)) end if end if end do IF ((ExcitMat3(1, 2) == 0) .and. (.not. tAllExcitFound)) then CALL FindNewSingDet(nI, nJ, OrbiIndex, OrbA, ExcitMat3, tParity) end if ENDSUBROUTINE GenSingleExcit SUBROUTINE GenDoubleExcit(nI, iLut, nJ, ExcitMat3, tParity, tAllExcitFound, tij_lt_ab_only) ! This generates one by one, all possible double excitations. ! This involves a way of ordering the electron pairs i,j and a,b so that given an i,j and a,b we can find the next. ! The overall symmetry must also be maintained - i.e. if i and j are alpha and beta, a and b must be alpha and beta ! or vice versa. USE SystemData, only: ElecPairs, tFixLz, iMaxLz USE GenRandSymExcitNUMod, only: PickElecPair use constants, only: bits_n_int INTEGER :: nI(NEl), Orbj, Orbi, Orba, Orbb, Syma, Symb INTEGER(KIND=n_int) :: iLut(0:NIfTot) INTEGER :: Elec1Ind, Elec2Ind, SymProduct, iSpn, Spinb, nJ(NEl), ExcitMat3(2, 2), SumMl INTEGER, SAVE :: ijInd, OrbaChosen, OrbbIndex, Spina, SymInd LOGICAL :: tDoubleExcitFound, tFirsta, tFirstb, tNewij, tNewa, tAllExcitFound, tParity, tij_lt_ab_only INTEGER :: Mla, Mlb, Indij tDoubleExcitFound = .false. tFirsta = .false. tFirstb = .false. Orbi = ExcitMat3(1, 1) Orbj = ExcitMat3(1, 2) Orba = ExcitMat3(2, 1) Orbb = ExcitMat3(2, 2) IF (Orbi == 0) THEN ijInd = 1 ! If Orbi, then we are choosing the first double. ! It is therefore also the first set of a and b for this electron pair i,j. tFirsta = .true. tFirstb = .true. end if lp: do while (.not. tDoubleExcitFound) ! Otherwise we use the previous ijInd and the saved indexes for a and b. ! This routine allows us to pick an electron pair i,j specified by the index ijInd. ! The i and j orbitals are then given by nI(Elec1Ind) and nI(Elec2Ind), and the symmetry product of the two is ! SymProduct and the spin iSpn. ! iSpn=2 for alpha beta pair, ispn=3 for alpha alpha pair and ispn=1 for beta beta pair. CALL PickElecPair(nI, Elec1Ind, Elec2Ind, SymProduct, iSpn, SumMl, ijInd) Indij = ((((nI(Elec2Ind) - 2) * (nI(Elec2Ind) - 1)) / 2) + nI(Elec1Ind)) tNewij = .false. ! This becomes true when we can no longer find an allowed orbital a for this ij pair and we need to move onto the next. do while ((.not. tNewij) .and. (.not. tDoubleExcitFound)) ! This loop runs through the allowed a orbitals ! until a double excitation is found. IF (tFirsta) THEN ! If this is the first double we are picking with this ij, we start with the alpha spin, unless i and j are both beta. ! There is no restriction on the symmetry for orbital a - although clearly the symmetry we pick determins b. ! write(stdout,*) "iSpn",iSpn IF (iSpn == 1) THEN !beta beta pair Spina = 2 !Want to start at the first beta orbital OrbaChosen = 1 else if (iSpn == 3) THEN !alpha alpha pair Spina = 1 OrbaChosen = 2 ELSE Spina = 2 OrbaChosen = 1 end if end if ! If it is not the first, we have stored the previous spina and orba index - need to start with these and see ! if any more double remain. Orba = OrbaChosen ! write(stdout,*) "Chosen index, orbital for a: ",OrbaChosen,Orba ! The orbital chosen must be unoccupied. This is just a test to make sure this is the case. do while ((BTEST(iLut((Orba - 1) / bits_n_int), MOD((Orba - 1), bits_n_int))) .or. (abs(SumMl - G1(Orba)%Ml) > iMaxLz)) !We also test that the summl value and ml of Orba is such that it is possible for orbb to conserve ml. !Will get into this loop if the orbital is occupied, or if the ml is such that no orbb is possible. ! If not, we move onto the next orbital. IF (iSpn /= 2) THEN !Increment by two, since we want to look at the same spin state. OrbaChosen = OrbaChosen + 2 ELSE !Increment by one, since we want to look at both alpha and beta spins. OrbaChosen = OrbaChosen + 1 IF (Spina == 2) THEN Spina = 1 ELSE Spina = 2 end if end if IF (OrbaChosen > nBasis) THEN !We have reached the end of all allowed symmetries for the a orbital, only taking !into account spin symmetry. Choose new ij pair now. tNewij = .true. EXIT end if ! Otherwise the new orbital a is the first unoccupied orbital of allowed symmetry etc. Orba = OrbaChosen ! write(stdout,*) "Chosen index, orbital for a: ",OrbaChosen,Orba end do ! If we have got to the end of the a orbitals, and need a new i,j pair, we increment ijInd and check ! this hasn't gone beyond the limits - bail out if we have. IF (tNewij) THEN ijInd = ijInd + 1 IF (ijInd > ElecPairs) THEN tDoubleExcitFound = .true. tAllExcitFound = .true. ! AllExcitFound true indicates there are no more symmetry allowed double excitations. end if EXIT end if tNewa = .false. !Find a b do while ((.not. tNewa) .and. (.not. tDoubleExcitFound)) ! We now have i,j,a and we just need to pick b. ! First find the spin of b. IF (iSpn == 1) THEN Spinb = 2 else if (iSpn == 3) THEN Spinb = 1 ELSE IF (Spina == 1) THEN Spinb = 2 ELSE Spinb = 1 end if end if ! Then find the symmetry of b. IF (tNoSymGenRandExcits) THEN Syma = 0 ELSE Syma = INT(G1(Orba)%Sym%S, 4) end if Symb = IEOR(Syma, SymProduct) ! Then find the ml of b. IF (tFixLz) THEN Mla = G1(Orba)%Ml Mlb = SumMl - Mla ELSE Mla = 0 Mlb = 0 end if ! If this is the first time we've picked an orbital b for these i,j and a, begin at the start of the symmetry block. ! Otherwise pick up where we left off last time. IF (tFirstb) THEN SymInd = ClassCountInd(Spinb, Symb, Mlb) OrbbIndex = SymLabelCounts2(1, SymInd) ELSE !Update new orbital b index OrbbIndex = OrbbIndex + 1 end if ! If the new b orbital is still within the limits, check it is unoccupied and move onto the next orbital if it is. IF (OrbbIndex > nBasis) THEN tNewa = .true. tFirsta = .false. end if IF (.not. tNewa) THEN IF (OrbbIndex > (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) THEN ! If we have already gone beyond the symmetry limits by choosing the next b orbital, pick a new a orbital. tNewa = .true. tFirsta = .false. ELSE Orbb = SymLabelList2(OrbbIndex) ! Checking the orbital b is unoccupied and > a. do while (((BTEST(iLut((Orbb - 1) / bits_n_int), MOD((Orbb - 1), bits_n_int))) .or. (Orbb <= Orba)) .or. & (tij_lt_ab_only .and. (((((Orbb - 2) * (Orbb - 1)) / 2) + Orba) < Indij))) !Orbital is occupied - try again OrbbIndex = OrbbIndex + 1 IF (OrbbIndex > (SymLabelCounts2(1, SymInd) + SymLabelCounts2(2, SymInd) - 1)) THEN !Reached end of symmetry block - need new a ! write(stdout,*) "Reached end of sym block",Orbb,Orba tNewa = .true. tFirsta = .false. EXIT end if ! write(stdout,*) "Cycling through orbitals: ",OrbbIndex,Orbb !Update new orbital b index Orbb = SymLabelList2(OrbbIndex) ! write(stdout,*) "Attempting again with orbital: ",Orbb end do end if end if ! If we are moving onto the next a orbital, check we don't also need a new ij pair. IF (tNewa) THEN IF (iSpn /= 2) THEN !Increment by two, since we want to look at the same spin state. OrbaChosen = OrbaChosen + 2 ! tFirsta=.false. ELSE !Increment by one, since we want to look at both alpha and beta spins. IF (Spina == 1) THEN Spina = 2 ELSE Spina = 1 end if OrbaChosen = OrbaChosen + 1 ! tFirsta=.false. end if tFirstb = .true. ! write(stdout,*) "New OrbaChosen: ",OrbaChosen IF (OrbaChosen > nBasis) THEN !We have reached the end of all allowed symmetries for the a orbital, only taking !into account spin symmetry. Choose new ij pair now. tNewij = .true. ijInd = ijInd + 1 ! write(stdout,*) "ijInd: ",ijInd IF (ijInd > ElecPairs) THEN tAllExcitFound = .true. tDoubleExcitFound = .false. EXIT lp end if end if ELSE !If we don't need a new a, we have found an excitation ij -> ab that is accepted. tDoubleExcitFound = .true. end if end do end do ! This is the loop for new ij pairs - if we are choosing a new ij we are automatically choosing a new a and b also. tFirsta = .true. tFirstb = .true. end do lp if (tDoubleExcitFound .and. (.not. tAllExcitFound)) then call make_double(nI, nJ, elec1ind, elec2ind, orbA, orbB, & ExcitMat3, tParity) end if ENDSUBROUTINE GenDoubleExcit !This routine creates the final determinant for a single excitation. SUBROUTINE FindNewSingDet(nI, nJ, Elec1Ind, OrbA, ExcitMat3, tParity) INTEGER :: nI(NEl), nJ(NEl), Elec1Ind, OrbA, ExcitMat3(2, 2) LOGICAL :: tParity !First construct ExcitMat3 ExcitMat3(1, 1) = Elec1Ind ExcitMat3(2, 1) = OrbA ExcitMat3(1, 2) = 0 ExcitMat3(2, 2) = 0 nJ(:) = nI(:) CALL FindExcitDet(ExcitMat3, nJ, 1, tParity) END SUBROUTINE FindNewSingDet !> @brief !> Return all configurations that are connected to nI as !> array of iluts (det_list(0:niftot, n_excits)). !> !> @details !> Triple excitations are not supported. !> !> @param[in] nI, The configuration from which to excite. !> @param[out] n_excits, The number of connected configurations. !> @param[out] det_list, The connected configurations in ilut format. !> (det_list(0:niftot, n_excits)) !> @param[in] ex_flag, The requested excitations. (1 = singles, 2 = doubles) !> If ommited all excitations will be generated. subroutine gen_excits(nI, n_excits, det_list, ex_flag) integer, intent(in) :: nI(nel) integer, intent(out) :: n_excits integer(n_int), intent(out), allocatable :: det_list(:, :) integer, optional, intent(in) :: ex_flag character(*), parameter :: this_routine = "gen_all_excits_default" integer :: n_singles, n_doubles, n_dets, ex(2, maxExcit), ex_flag_ integer :: nJ(nel) logical :: tpar, found_all integer(n_int) :: ilut(0:niftot) integer, parameter :: arbitrary_number = 42 ! can be neither 1 or 2 n_excits = -1 call EncodeBitDet(nI, ilut) if (present(ex_flag)) then ASSERT(present(ex_flag) .implies. any(ex_flag == [1, 2])) end if ! If it is set to neither 1 nor 2, all excitations ! are generated def_default(ex_flag_, ex_flag, arbitrary_number) ! for reference in the "normal" case it looks like that: call CountExcitations3(nI, ex_flag_, n_singles, n_doubles) n_excits = n_singles + n_doubles allocate(det_list(0:niftot, n_excits)) n_dets = 0 ex = 0 call GenExcitations3(nI, ilut, nJ, ex_flag_, ex, tpar, found_all, & .false.) do while (.not. found_all) n_dets = n_dets + 1 call EncodeBitDet(nJ, det_list(:, n_dets)) call GenExcitations3(nI, ilut, nJ, ex_flag_, ex, tpar, & found_all, .false.) end do if (n_dets /= n_excits) then write(stdout, *) "expected number of excitations: ", n_excits write(stdout, *) "actual calculated ones: ", n_dets call stop_all(this_routine, "Incorrect number of excitations found") end if ! Sort the dets, so they are easy to find by binary searching call sort(det_list, ilut_lt, ilut_gt) end subroutine gen_excits !> @brief !> Return all configurations that are connected to nI !> as array of iluts (det_list(0:niftot, n_excits)). subroutine gen_all_excits(nI, n_excits, det_list) integer, intent(in) :: nI(nel) integer, intent(out) :: n_excits integer(n_int), intent(out), allocatable :: det_list(:, :) call gen_excits(nI, n_excits, det_list) end subroutine gen_all_excits END MODULE SymExcit3