calcDiagMatEleGuga_nI Function

    pure function calcDiagMatEleGuga_nI(nI) result(hel_ret)
        ! calculates the diagonal Hamiltonian matrix element when a CSF in
        ! nI(nEl) form is provided and returns hElement of type hElement_t
        integer, intent(in) :: nI(nEl)
        HElement_t(dp) :: hel_ret

        ! have to loop over the number of spatial orbitals i , and within
        ! loop again over orbitals j > i, s indicates spatial orbitals
        integer :: iOrb, jOrb, inc1, inc2, sOrb, pOrb
        real(dp) :: nOcc1, nOcc2

        hel_ret = ECore

        iOrb = 1
        ! loop over nI spin orbital entries: good thing is unoccupied orbitals
        ! do not  contribute to the single matrix element part.
        do while (iOrb <= nEl)
            ! spatial orbital index needed for get_umat_el access
            sOrb = (nI(iOrb) + 1) / 2
            ! have to check if orbital is singly or doubly occupied.
            if (isDouble(nI, iOrb)) then ! double has two part. int. contribution
                nOcc1 = 2.0_dp
                hel_ret = hel_ret + nOcc1 * GetTMatEl(nI(iOrb), nI(iOrb)) + &
                          get_umat_el(sOrb, sOrb, sOrb, sOrb)

                ! correctly count through spin orbitals if its a double occ.
                inc1 = 2

            else ! single occupation
                nOcc1 = 1.0_dp
                hel_ret = hel_ret + nOcc1 * GetTMatEl(nI(iOrb), nI(iOrb))
                inc1 = 1

            end if

            ! second loop:
            jOrb = iOrb + inc1
            do while (jOrb <= nEl)
                pOrb = (nI(jOrb) + 1) / 2
                ! again check for double occupancies
                if (isDouble(nI, jOrb)) then
                    nOcc2 = 2.0_dp
                    inc2 = 2

                else
                    nOcc2 = 1.0_dp
                    inc2 = 1

                end if
                ! standard two particle contribution
                if (.not. (t_tJ_model .or. t_heisenberg_model)) then
                    hel_ret = hel_ret + nOcc1 * nOcc2 * ( &
                              get_umat_el(sOrb, pOrb, sOrb, pOrb) - &
                              get_umat_el(sOrb, pOrb, pOrb, sOrb) / 2.0_dp)
                end if

                ! calculate exchange integral part, involving Shavitt
                ! rules for matrix elements, only contributes if both
                ! stepvectors are 1 or 2, still to be implemented..
                if ((nOcc1.isclose.1.0_dp) .and. (nOcc2.isclose.1.0_dp)) then
                    hel_ret = hel_ret - get_umat_el(sOrb, pOrb, pOrb, sOrb) * &
                              calcDiagExchangeGUGA_nI(iOrb, jOrb, nI) / 2.0_dp

                    if (t_tJ_model) then
                        hel_ret = hel_ret + get_umat_el(sOrb, pOrb, pOrb, sOrb) / 2.0
                    end if
                end if

                ! increment the counters
                jOrb = jOrb + inc2

            end do
            iOrb = iOrb + inc1
        end do

    end function calcDiagMatEleGUGA_nI