Par2vSum Subroutine

    Subroutine Par2vSum(nI)
        !=  A parallel version of the 2-vertex sum.
        !=
        !=   In:
        !=     nI(nEl)        The root determinant of the 2v sum.
        !=
        !=  This is not quite stable yet.
        !=
        !=  Issues:
        !=     * Some problems remain with how electrons are distributed over processors.
        !=     * Doesn't work for CPMD calculations.
        use constants, only: dp
        use SystemData, only: nEl, Beta
        Use Determinants, only: get_helement
        use neci_intfce
        IMPLICIT NONE
        Integer nI(nEl)
        integer iMinElec, iMaxElec
        integer i
        integer store(6)
        integer ic, exlen(1), iC0
        integer, pointer :: Ex(:)
        integer nJ(nEl)
        HElement_t(dp) dU
        real(dp) dE1, dE2
        HElement_t(dp) dEw, dw, dEwtot, dwtot, dTots(2), dTots2(2)
        iC0 = 0
        i = iProcIndex + 1
        write(stdout, *) "Proc ", i, "/", nProcessors
        call GetProcElectrons(iProcIndex + 1, iMinElec, iMaxElec)
        write(stdout, *) "Electrons ", iMinElec, " TO ", iMaxElec

!  The root's energy
        dE1 = get_helement(nI, nI, 0)

!  Initialize.  If we're the first processor then we add in the 1-vertex graph.
        if (iProcIndex == 0) THEN
            dEwTot = dE1
            dwTot = (1.0_dp)
        ELSE
            dEwTot = 0.0_dp
            dwTot = (0.0_dp)
        end if

! Now enumerate all 2v graphs
!.. Setup the spin excit generator
        STORE(1) = 0
!  IC is the excitation level (relative to the reverence det).
        CALL GENSYMEXCITIT3Par(NI, .TRUE., EXLEN, nJ, IC, STORE, 3, iMinElec, iMaxElec)
        allocate(Ex(exLen(1)))
        EX(1) = 0
        CALL GENSYMEXCITIT3Par(NI, .TRUE., EX, nJ, IC, STORE, 3, iMinElec, iMaxElec)

!  Generate the first excitation
        CALL GENSYMEXCITIT3Par(NI, .False., EX, nJ, IC, STORE, 3, iMinElec, iMaxElec)
        i = 0
!NJ(1) is zero when there are no more excitations.
        DO WHILE (NJ(1) /= 0)
            i = i + 1
            dU = get_helement(nI, nJ, IC) ! Whilst we know IC, we don't know tSign!
            dE2 = get_helement(nJ, nJ, 0)
            call Get2vWeightEnergy(dE1, dE2, dU, Beta, dw, dEw)
            dEwTot = dEwTot + dEw
            dwTot = dwTot + dw
            CALL GENSYMEXCITIT3Par(NI, .False., EX, nJ, IC, STORE, 3, iMinElec, iMaxElec)
        end do
        write(stdout, *) I
        write(stdout, *) dEwTot, dwTot, dEwTot / dwTot
        dTots(1) = dwTot
        dTots(2) = dEwTot
        Call MPISumAll(dTots, 2, dTots2)
        write(stdout, *) dTots2(2), dTots2(1), dTots2(2) / dTots2(1)
        deallocate(Ex)
    End Subroutine Par2vSum