FindNewDet Subroutine

public subroutine FindNewDet(session, nJ, tParity)

Arguments

Type IntentOptional Attributes Name
type(ExcitGenSessionType), intent(in) :: session
integer, intent(out) :: nJ(nEl)
logical, intent(out) :: tParity

Contents

Source Code


Source Code

    subroutine FindNewDet(session, nJ, tParity)
        implicit none
        type(ExcitGenSessionType), intent(in) :: session
        integer, intent(out) :: nJ(nEl)
        logical, intent(out) :: tParity
        integer :: nJtmp(nEL)
        integer :: minOrbBound, maxOrbBound
        integer :: holes(nBasis - nEl)
        integer :: i, j, switchedElecs
        integer :: elecIdx, elecOrb, holeOrb
        nJ = session%nI
        holes = session%holes

        tParity = .false. ! even parity

        do i = 1, session%rank
            ! loop over electron hole pairs
            ! spin orb arrays in session do the same job as an excitmat

            elecOrb = session%elecSpinOrbs(i)
            holeOrb = session%holeSpinOrbs(i)

            minOrbBound = min(elecOrb, holeOrb)
            maxOrbBound = max(elecOrb, holeOrb)
            switchedElecs = 0
            do j = 1, nEl
                ! the number of electrons inbetween determines parity
                ! and also determines insertion position
                if (nJ(j) > minOrbBound .and. nJ(j) < maxOrbBound) then
                    switchedElecs = switchedElecs + 1
                    tParity = .not. tParity
                end if

            end do

            ! Suppose we have the following excitation
            !
            !   1  2  3  4  5  6  7  8   (spinorbs)
            !
            ! [ 1, 0, 1, 1, 1, 0, 0, 0 ]
            !              |
            !              V
            ! [ 1, 1, 0, 0, 1, 0, 1, 0 ]
            !
            ! this is a rank 2 excitation via:
            ! elecIndices = (2, 3)
            ! holeIndices = (3, 1)
            !
            ! which gives:
            ! elecSpinOrbs = (3, 4)
            ! holeSpinOrbs = (7, 2)
            !
            ! the first excitation is 3 -> 7
            !
            ! nJ was originally:
            ! ( 1, 3, 4, 5 )
            ! first iteration we should get:
            ! ( 1, 4, 5, 7 )
            ! this routine must yield:
            ! ( 1, 2, 5, 7 )
            !
            ! the loop above has given us the number of electrons
            ! interchanged in the excitation
            !
            ! we have all the information we need to manipulate nJ
            ! such that it is still ordered after the first excitation
            !
            ! for a hole orbital greater than the electron orbital
            ! the insertion point is at elecIdx + switchedElecs
            !
            ! 2 switched elecs
            ! use a temporary array nJtmp
            !
            ! business as usual for the part of the array before the moving electron
            ! nJtmp(1:elecIdx-1) = nJ(1:elecIdx-1)
            !    ( 1, 0, 0, 0 )
            !
            ! move down the elements between origin and destination indices
            ! nJtmp(elecIdx:elecIdx+switchedElecs-1) = nJ(elecIdx+1:elecIdx+switchedElecs)
            !    ( 1, 4, 5, 0 )
            !
            ! insert the electron at its destination orbital
            ! nJtmp(elecIdx+switchedElecs) = holeOrb
            !    ( 1, 4, 5, 7 )
            !
            ! finally, simply copy the rest of the array
            ! nJtmp(elecIdx+switchedElecs+1:nEl) = nJ(elecIdx+switchedElecs+1:nEl)
            !    ( 1, 4, 5, 7 )
            !
            !
            ! the second excitation is 4 -> 2
            !
            ! likewise, for an electron orbital greater than the hole
            ! orbital the insertion point is at elecIdx - switchedElecs
            !
            ! but there are 0 switched elements, so simply do
            ! nJ(elecIdx) = holeOrb

            ! we can't just take elecIdx from session%elecIndices, because of the re-ordering
            ! that occurs in this loop, so first we have to find it

            do j = 1, nEl
                if (nJ(j) == elecOrb) then
                    elecIdx = j
                end if
            end do

            ! then construct the arrays
            if (switchedElecs == 0) then
                ! this is easy
                nJ(elecIdx) = holeOrb
            else
                if (holeOrb - elecOrb > 0) then
                    nJtmp(1:elecIdx - 1) = nJ(1:elecIdx - 1)
                    nJtmp(elecIdx:elecIdx + switchedElecs - 1) = nJ(elecIdx + 1:elecIdx + switchedElecs)
                    nJtmp(elecIdx + switchedElecs) = holeOrb
                    nJtmp(elecIdx + switchedElecs + 1:nEl) = nJ(elecIdx + switchedElecs + 1:nEl)
                else
                    ! same as above, but with elecIdx -> elecIdx - switchedElecs
                    nJtmp(1:elecIdx - switchedElecs - 1) = nJ(1:elecIdx - switchedElecs - 1)
                    nJtmp(elecIdx - switchedElecs:elecIdx - 1) = nJ(elecIdx - switchedElecs + 1:elecIdx)
                    nJtmp(elecIdx) = holeOrb
                    nJtmp(elecIdx + 1:nEl) = nJ(elecIdx + 1:nEl)
                end if
                ! copy back into target det array
                nJ(1:nEl) = nJtmp
            end if
        end do
    end subroutine FindNewDet