init_rdm_estimates_t Subroutine

    subroutine init_rdm_estimates_t(est, nrdms_standard, nrdms_transition, open_output_file, &
                                    filename)

        ! Initialise an rdm_estimates_t object. Allocate arrays to be large
        ! enough to hold estimates for nrdms_srandard+nrdms_transition RDMs.

        ! Also, if open_output_file is true, and if this is the processor with
        ! label 0, then open an RDMEstimates file, and write the file's header.

        use CalcData, only: tEN2
        use LoggingData, only: tCalcPropEst, iNumPropToEst
        use Parallel_neci, only: iProcIndex
        use rdm_data, only: rdm_estimates_t
        use util_mod, only: get_free_unit

        type(rdm_estimates_t), intent(out) :: est
        integer, intent(in) :: nrdms_standard, nrdms_transition
        logical, intent(in) :: open_output_file
        character(*), intent(in), optional :: filename

        integer :: nrdms, ierr
        character(255) :: rdm_filename

        nrdms = nrdms_standard + nrdms_transition

        ! Store the number of RDMs.
        est%nrdms = nrdms
        est%nrdms_standard = nrdms_standard
        est%nrdms_transition = nrdms_transition

        ! Estimates over the entire RDM sampling period.
        allocate(est%trace(nrdms), stat=ierr)
        allocate(est%norm(nrdms), stat=ierr)
        allocate(est%energy_1_num(nrdms), stat=ierr)
        allocate(est%energy_2_num(nrdms), stat=ierr)
        allocate(est%energy_num(nrdms), stat=ierr)
        if (.not. tGUGA) allocate(est%spin_num(nrdms), stat=ierr)
        if (tCalcPropEst) allocate(est%property(iNumPropToEst, nrdms), stat=ierr)
        if (tEN2) allocate(est%energy_pert(nrdms_standard), stat=ierr)

        ! "Instantaneous" estimates over the previous sampling block.
        allocate(est%trace_inst(nrdms), stat=ierr)
        allocate(est%norm_inst(nrdms), stat=ierr)
        allocate(est%energy_1_num_inst(nrdms), stat=ierr)
        allocate(est%energy_2_num_inst(nrdms), stat=ierr)
        allocate(est%energy_num_inst(nrdms), stat=ierr)
        if (.not. tGUGA) allocate(est%spin_num_inst(nrdms), stat=ierr)
        if (tCalcPropEst) allocate(est%property_inst(iNumPropToEst, nrdms), stat=ierr)
        if (tEN2) allocate(est%energy_pert_inst(nrdms_standard), stat=ierr)

        ! Hermiticity errors, for the final RDMs.
        allocate(est%max_error_herm(nrdms), stat=ierr)
        allocate(est%sum_error_herm(nrdms), stat=ierr)

        est%trace = 0.0_dp
        est%norm = 0.0_dp
        est%energy_1_num = 0.0_dp
        est%energy_2_num = 0.0_dp
        est%energy_num = 0.0_dp
        if (.not. tGUGA) est%spin_num = 0.0_dp
        if (tCalcPropEst) est%property = 0.0_dp
        if (tEN2) est%energy_pert = 0.0_dp

        est%trace_inst = 0.0_dp
        est%norm_inst = 0.0_dp
        est%energy_1_num_inst = 0.0_dp
        est%energy_2_num_inst = 0.0_dp
        est%energy_num_inst = 0.0_dp
        if (.not. tGUGA) est%spin_num_inst = 0.0_dp
        if (tCalcPropEst) est%property_inst = 0.0_dp
        if (tEN2) est%energy_pert_inst = 0.0_dp

        est%max_error_herm = 0.0_dp
        est%sum_error_herm = 0.0_dp

        ! If appropriate, create a new RDMEstimates file.
        if (iProcIndex == 0 .and. open_output_file) then
            est%write_unit = get_free_unit()
            if (present(filename)) then
                rdm_filename = filename
            else
                rdm_filename = "RDMEstimates"
            end if
            call write_rdm_est_file_header(est%write_unit, nrdms_standard, nrdms_transition, &
                                           rdm_filename)
        else
            ! If we don't have a file open with this unit, set it to something
            ! unique, so we can easily check, and which will cause an obvious
            ! error if we don't.
            est%write_unit = huge(est%write_unit)
        end if

    end subroutine init_rdm_estimates_t