CalcDoCalc Subroutine

    subroutine CalcDoCalc(kp)
        use SystemData, only: Alat, Arr, Brr, Beta, ECore, G1, LMS, LMS2, nBasis, NMSH, nBasisMax
        use SystemData, only: SymRestrict, tParity, tSpn, ALat, Beta, tMolpro, tMolproMimic
        use SystemData, only: Symmetry, SymmetrySize, SymmetrySizeB, BasisFN, BasisFNSize, BasisFNSizeB, nEl
        Use DetCalcData, only: nDet, nEval, nmrks, w
        USE FciMCParMod, only: FciMCPar
        use RPA_Mod, only: RunRPA_QBA
        use DetCalc, only: CK, DetInv, tEnergy, tRead
        use IntegralsData, only: FCK, NMAX, UMat, FCK
        use IntegralsData, only: HFEDelta, HFMix, nTay
        Use LoggingData, only: iLogging
        use Parallel_Calc
        use util_mod, only: get_free_unit, NECI_ICOPY
        use sym_mod
        use davidson_neci, only: DavidsonCalcType, DestroyDavidsonCalc
        use davidson_neci, only: davidson_direct_ci_init, davidson_direct_ci_end, perform_davidson
        use hamiltonian_linalg, only: direct_ci_type
        use kp_fciqmc, only: perform_kp_fciqmc, perform_subspace_fciqmc
        use kp_fciqmc_data_mod, only: tExcitedStateKP
        use kp_fciqmc_procs, only: kp_fciqmc_data
        use util_mod, only: int_fmt

        real(dp) :: EN, WeightDum, EnerDum
        real(dp), allocatable :: final_energy(:)
        integer :: iSeed, iunit, i
        type(kp_fciqmc_data), intent(inout) :: kp
        character(*), parameter :: this_routine = 'CalcDoCalc'
        type(DavidsonCalcType) :: davidsonCalc

        iSeed = 7

        IF(tMP2Standalone) then
            call ParMP2(FDet)
            ! Parallal 2v sum currently for testing only.
            !          call Par2vSum(FDet)
        ELSE IF(tDavidson) then
            davidsonCalc = davidson_direct_ci_init()
            if(t_non_hermitian_2_body) then
                call stop_all(this_routine, &
                              "perform_davidson not adapted for non-hermitian Hamiltonians!")
            end if
            if(tGUGA) then
                call stop_all(this_routine, &
                              "perform_davidson not adapted for GUGA yet")
            end if
            call perform_davidson(davidsonCalc, direct_ci_type, .true.)
            call davidson_direct_ci_end(davidsonCalc)
            call DestroyDavidsonCalc(davidsonCalc)
        else if(allocated(pgen_unit_test_spec)) then
            call batch_run_excit_gen_tester(pgen_unit_test_spec)

        ELSE IF(NPATHS /= 0 .OR. DETINV > 0) THEN
            !Old and obsiolecte
            !             IF(TRHOIJND) THEN
            !C.. We're calculating the RHOs for interest's sake, and writing them,
            !C.. but not keeping them in memory
            !                  write(stdout,*) "Calculating RHOS..."
            !                  write(stdout,*) "Using approx NTAY=",NTAY
            !                  CALL CALCRHOSD(NMRKS,BETA,I_P,I_HMAX,I_VMAX,NEL,NDET,        &
            !     &               NBASISMAX,G1,nBasis,BRR,NMSH,FCK,NMAX,ALAT,UMAT,             &
            !     &               NTAY,RHOEPS,NWHTAY,ECORE)
            !             end if

            if(tFCIMC) then
                call FciMCPar(final_energy)
                if((.not. tMolpro) .and. (.not. tMolproMimic)) then
                    if(allocated(final_energy)) then
                        do i = 1, size(final_energy)
                            write(stdout, '(1X,"Final energy estimate for state",1X,'//int_fmt(i)//',":",g25.14)') &
                                i, final_energy(i)
                        end do
                    end if
                end if
            else if(tRPA_QBA) then
                call RunRPA_QBA(WeightDum, EnerDum)
                write(stdout, *) "Summed approx E(Beta)=", EnerDum
            else if(tKP_FCIQMC) then
                if(tExcitedStateKP) then
                    call perform_subspace_fciqmc(kp)
                else
                    call perform_kp_fciqmc(kp)
                end if
            else if(tRPA_QBA) then
                call RunRPA_QBA(WeightDum, EnerDum)
                write(stdout, *) "Summed approx E(Beta)=", EnerDum
            else if(tKP_FCIQMC) then
                if(tExcitedStateKP) then
                    call perform_subspace_fciqmc(kp)
                else
                    call perform_kp_fciqmc(kp)
                end if
            else if(t_real_time_fciqmc) then
                call perform_real_time_fciqmc()
            end if
            IF(TMONTE .and. .not. tMP2Standalone) THEN
!             DBRAT=0.01
!             DBETA=DBRAT*BETA
                write(stdout, *) "I_HMAX:", I_HMAX
                write(stdout, *) "Calculating MC Energy..."
                CALL neci_flush(stdout)
                IF(NTAY(1) > 0) THEN
                    write(stdout, *) "Using approx RHOs generated on the fly, NTAY=", NTAY(1)
!C.. NMAX is now ARR
                    call stop_all(this_routine, "DMONTECARLO2 is now non-functional.")
                else if(NTAY(1) == 0) THEN
                    IF(TENERGY) THEN
                        write(stdout, *) "Using exact RHOs generated on the fly"
!C.. NTAY=0 signifying we're going to calculate the RHO values when we
!C.. need them from the list of eigenvalues.
!C.. Hide NMSH=NEVAL
!C..         FCK=W
!C..         ZIA=CK
!C..         UMAT=NDET
!C..         ALAT=NMRKS
!C..         NMAX=ARR
                        call stop_all(this_routine, "DMONTECARLO2 is now non-functional.")
!                   EN=DMONTECARLO2(MCDET,I_P,BETA,DBETA,I_HMAX,I_VMAX,IMCSTEPS,             &
!     &                G1,NEL,NBASISMAX,nBasis,BRR,IEQSTEPS,                                 &
!     &                NEVAL,W,CK,ARR,NMRKS,NDET,NTAY,RHOEPS,NWHTAY,ILOGGING,ECORE,BETAEQ)
                    ELSE
                        call stop_all(this_routine, "TENERGY not set, but NTAY=0")
                    end if
                end if
                write(stdout, *) "MC Energy:", EN
!CC           write(12,*) DBRAT,EN
            end if
        end if
!C.. /AJWT
    End Subroutine CalcDoCalc