PrintOrbOccs Subroutine

    subroutine PrintOrbOccs(OrbOccs)

        ! This routine takes whatever orbital basis we're using and is called
        ! at the end of a spawn to find the contribution of each orbital to the
        ! final wavefunction. This is done by histogramming the determinant
        ! populations, and then running over these adding the  coefficients of
        ! each determinant to the orbitals occupied. This is essentially
        ! < Psi | a_p+ a_p | Psi > - the diagonal terms of the one electron
        ! reduced density matrix.

        real(dp) :: Norm, OrbOccs(nBasis), AllOrbOccs(nBasis)
        integer :: i, iunit
        logical :: tWarning

        AllOrbOccs = 0.0_dp

        call MPIReduce(OrbOccs, MPI_SUM, AllOrbOccs)

        ! Want to normalise the orbital contributions for convenience.
        tWarning = .false.
        if (iProcIndex == 0) then
            Norm = 0.0_dp
            do i = 1, nBasis
                Norm = Norm + AllOrbOccs(i)
                if ((AllOrbOccs(i) < 0) .or. (Norm < 0)) then
                    write(stdout, *) 'WARNING: Integer overflow when calculating the orbital occupations.'
                    tWarning = .true.
                end if
            end do
            if (Norm > 1.0e-8_dp) then
                do i = 1, nBasis
                    AllOrbOccs(i) = AllOrbOccs(i) / Norm
                end do
            end if

            iunit = get_free_unit()
            open(iunit, file='ORBOCCUPATIONS', status='UNKNOWN')
            write(iunit, '(A15,A30)') '# Orbital no.', 'Normalised occupation'
            if (tWarning) write(iunit, *) 'WARNING: integer OVERFLOW OCCURRED WHEN CALCULATING THESE OCCUPATIONS'
            do i = 1, nBasis
                write(iunit, '(I15,F30.10)') i, AllOrbOccs(i)
            end do
            close(iunit)
        end if

    end subroutine PrintOrbOccs