FindNatOrbs Subroutine

    subroutine FindNatOrbs()

        ! Fed into this routine will be the wavefunction and its amplitudes
        ! within the given excitation level.

        ! First need to set up the orbital labels and symmetries etc. This is
        ! done slightly differently for spin and spatial and whether or not we
        ! are truncating the virtual space when writing out the final ROFCIDUMP
        ! file.

        ! Allocate the matrix used to find the natural orbitals.

        integer :: ierr
        character(len=*), parameter :: t_r = 'FindNatOrbs'

        allocate(NatOrbMat(NoOrbs, NoOrbs), stat=ierr)
        call LogMemAlloc('NatOrbMat', NoOrbs**2, 8, t_r, NatOrbMatTag, ierr)
        if (ierr /= 0) call stop_all(t_r, "Mem allocation for NatOrbMat failed.")
        NatOrbMat(:, :) = 0.0_dp

        allocate(Evalues(NoOrbs), stat=ierr)
        call LogMemAlloc('Evalues', NoOrbs, 8, t_r, EvaluesTag, ierr)
        if (ierr /= 0) call stop_all(t_r, "Mem allocation for Evalues failed.")
        Evalues(:) = 0.0_dp

        ! First need to fill the relevant matrix for calculating the type of
        ! natural orbitals we want.
        if (tFindCINatOrbs) then

            ! For the CISD, CISDT etc natural orbitals, the relevant matrix is
            ! the one electron reduced density matrix from the previous
            ! spawning calculation (trucated at a certain excitation).
            call FillOneRDM()

        else if (tUseMP2VarDenMat) then

            ! For the MP2 natural orbitals, the natural orbital matrix is the
            ! MP2 variational density matrix.
            call FillMP2VDM()

        end if

        ! Then need to diagonalise this, maintaining the various symmetries
        ! (spin and spatial).
        call DiagNatOrbMat()

        ! We then need to put the resulting eigenvectors back into the ordering
        ! we want, and copy these over to CoeffT1.
        call OrderCoeffT1()

    end subroutine FindNatOrbs