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), walker_norm
        integer :: i, iunit
        logical :: tWarning
        integer(n_int) , allocatable :: GlobalLargestWalkers(:, :)
        integer :: this_run
        real(dp), dimension(lenof_sign) :: currsign
        real(dp) :: currsignreal

        allocate(GlobalLargestWalkers(0:NIfTot,nBasis), source=0_n_int)

        AllOrbOccs = 0.0_dp

        call MPIReduce(OrbOccs, MPI_SUM, AllOrbOccs)
        this_run = GLOBAL_RUN
        call global_most_populated_states(nBasis, this_run, GlobalLargestWalkers, norm=walker_norm)

        ! Want to normalise the orbital contributions for convenience.
        tWarning = .false.
        if (iProcIndex == 0) then
            iunit = get_free_unit()
            ! normalise to nel
            Norm = 0.0_dp
            do i = 1, nBasis
                Norm = Norm + AllOrbOccs(i)

                if ((AllOrbOccs(i) < 0) .or. (Norm < 0)) then
                    write(stderr, *) '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) = nEl * AllOrbOccs(i) / Norm
                end do
            end if

            open(iunit, file='ORBOCCUPATIONS', status='UNKNOWN')
            write(iunit, '(A15,A30,A30,A30)') '# Orbital no.', 'occupation', 'Normalised CI coeff.', 'excitation'
            if (tWarning) write(iunit, *) 'WARNING: integer OVERFLOW OCCURRED WHEN CALCULATING THESE OCCUPATIONS'
            do i = 1, nBasis
                call extract_sign(GlobalLargestWalkers(:,i), currsign)
                currsignreal = core_space_weight(currsign, this_run)
                write(iunit, '(I15,F30.10,E30.12)', advance='no') i, AllOrbOccs(i), currsignreal/walker_norm
                write(iunit, '(A4)', advance='no') '    '
                call WriteDetBit(iunit,GlobalLargestWalkers(:,i),.true.)
            end do
            close(iunit)
        end if

    end subroutine PrintOrbOccs