Get2vWeightEnergy Subroutine

    subroutine Get2vWeightEnergy(dE1, dE2, dU, dBeta, dw, dEt)
        != Get the two-vertex contribution from a graph containing the reference
        != determinant, 0, and a connected determinant, i, given its matrix elements.
        != Weights are divided by exp(-beta E1)

        != Each two vertex graph is represented by a 2x2 (Hermitian) matrix
        !=  | dE1   dU  |
        !=  | dU*   dE2 |
        != The eigenvalues are [ (dE1+dE2) +- \Sqrt( (dE1+dE2) - 4dE1*dE2 + 4|dU|^2 ) ]/2.
        != where:
        !=   dE1=<D_0|H|D_0>
        !=   dE2=<D_i|H|D_i>
        !=   dU =<D_0|H|D_i>
        != Denoting the eigenvalues as dEp and dEm, the normalised eigenvectors are:
        !=   \frac{1}{\Sqrt{ dU^2 + (dE1-dEp)^2 } ( U, dEp-dE1 )
        !=   \frac{1}{\Sqrt{ dU^2 + (dE1-dEm)^2 } ( U, dEm-dE1 )

        != In:
        !=    dE1    <D_0|H|D_0>
        !=    dE2    <D_i|H|D_i>
        !=    dU     <D_0|H|D_i>
        !=    dBeta  beta
        != Out:
        !=    dw     weight of the graph, w_i[G]
        !=    dEt    weighted contribution of the graph, w_i[G] \tilde{E}_i[G]
        use constants, only: dp
        implicit none
        real(dp) dE1, dE2
        HElement_t(dp) dU, dEt, dw
        real(dp) dBeta
        HElement_t(dp) dEp, dEm, dD, dEx, dD2, dTmp
        if (near_zero(0.0_dp)) then
            ! Determinants are not connected.
            ! => zero contribution.
            dw = 0.0_dp
            dEt = 0.0_dp
            return
        end if

!  Calculate eigenvalues.
!   write(stdout,*) dE1,dE2,dU

        dD = ((dE1 + dE2)**2 - 4 * (dE1 * dE2) + 4 * abs(dU)**2)

!   write(stdout,*) dD

        dD = sqrt(dD) / 2
        dEp = (dE1 + dE2) / 2
        dEm = dEp - dD
        dEp = dEp + dD

!   write(stdout,*) dD,dEp,dEm

!  The normalized first coefficient of an eigenvector is U/sqrt((dE1-dEpm)**2+U**2)
        dD = 1 / sqrt((dE1 - dEp)**2 + abs(dU)**2) ! normalisation factor
        dD2 = dD * (dEp - dE1)               ! second coefficient
        dD = dD * dU                           ! first coefficient

!   write(stdout,*) dD,dD2

!dD is the eigenvector component
        dEx = exp(-dBeta * (dEp - dE1))
        dw = abs(dD)**2 * dEx
        dEt = dE1 * abs(dD)**2 * dEx
#ifdef CMPLX_
        dEt = dEt + dU * dD * conjg(dD2) * dEx
#else
        dEt = dEt + dU * dD * (dD2) * dEx
#endif

!   write(stdout,*) dEx,dw,dEt

!   write(stdout,*) dEp,dD,dD2,dw,dEx,dBeta
!  This can be numerically unstable when dE1 is v close to dEm:
!      dD=1/sqrt((dE1-dEm)**2+dU**2)
!      dD2=dD*(dEm-dE1)
!      dD=dD*dU
!  Instead we just swap dD2 and dD around
        dTmp = dD
#ifdef CMPLX_
        dD = conjg(dD2)
        dD2 = -conjg(dTmp)
#else
        dD = (dD2)
        dD2 = -(dTmp)
#endif
        dEx = exp(-dBeta * (dEm - dE1))
        dw = dw + (abs(dD)**2 * dEx)
!   write(stdout,*) dEm,dD,dD2,dw,dEx
        dEt = dEt + (dE1 * abs(dD)**2 * dEx)
#ifdef CMPLX_
        dEt = dEt + dU * dD * conjg(dD2) * dEx
#else
        dEt = dEt + dU * dD * (dD2) * dEx
#endif

!  Now, we already have calculated the effect of the one-vertex graph, so
!  we need to remove the double-counting of this from the 2-vertex graph.
!  For the one-vertex graph containing just i:
!     w_i[i] = exp(-\beta E_i)
!     w_i[i] \tilde{E}_i[i] = exp(-\beta E_i) E_i
!  As we factor out exp(-\beta E_i), this becomes just:
        dEt = dEt - (dE1)
        dw = dw - (1)
        write(stdout, *) 'wE,E', dEt, dw

    end subroutine Get2vWeightEnergy