DoDetCalc Subroutine

public subroutine DoDetCalc()

Arguments

None

Contents

Source Code


Source Code

    Subroutine DoDetCalc
        Use global_utilities
        use util_mod, only: get_free_unit
        use SystemData, only: Alat, arr, brr, boa, box, coa, ecore, g1, Beta
        use SystemData, only: t_new_real_space_hubbard
        use SystemData, only: nBasis, nBasisMax, nEl, nMsh, LzTot, TSPN, LMS
        use IntegralsData, only: FCK, NMAX, UMat
        Use LoggingData, only: iLogging, tLogDets, tCalcVariationalEnergy
        use Parallel_neci, only: iProcIndex
        use DetBitops, only: DetBitEQ, EncodeBitDet, FindBitExcitLevel
        use bit_rep_data, only: nifd, NIfTot, NIfD
        use bit_reps, only: decode_bit_det
        use sym_mod
        use HElem
        use MemoryManager, only: TagIntType
        use hist_data, only: tHistSpawn
        use CalcData, only: tFCIDavidson
        !use davidson_neci, only: DavidsonCalcType, DestroyDavidsonCalc
        !use davidson_neci, only: perform_davidson
        !use sparse_arrays, only: calculate_sparse_hamiltonian, deallocate_sparse_ham, &
        !sparse_ham, hamil_diag, HDiagTag, SparseHamilTags
        !use hamiltonian_linalg, only: sparse_hamil_type
        !use constants, only: size_n_int

        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, tCalcHFIndex
        use FCIMCData, only: davidson_ras, davidson_classes
        use ras, only: generate_entire_ras_space
        use real_space_hubbard, only: init_real_space_hubbard

        real(dp), ALLOCATABLE :: TKE(:), A(:, :), V(:), AM(:), BM(:), T(:), WT(:), SCR(:), WH(:), V2(:, :), FCIGS(:)
        HElement_t(dp), ALLOCATABLE :: WORK(:), WORK2(:)
        INTEGER, ALLOCATABLE :: INDEX(:), ISCR(:), Temp(:)
        integer(TagIntType) :: TKETag = 0, ATag = 0, VTag = 0, AMTag = 0, BMTag = 0, TTag = 0
        INTEGER(TagIntType) :: WTTag = 0, SCRTag = 0, ISCRTag = 0, INDEXTag = 0, WHTag = 0, Work2Tag = 0, V2Tag = 0, WorkTag = 0
        integer :: ierr, Lz
        character(25), parameter :: this_routine = 'DoDetCalc'
        real(dp) EXEN, GSEN
        Type(BasisFn) ISym, IHFSYM
        INTEGER GC, I, ICMAX
        INTEGER IN, IND, INDZ
        INTEGER NBLOCK!,OpenOrbs,OpenOrbsSym,Ex(2,NEl)
        INTEGER nKry1
        INTEGER(KIND=n_int) :: ilut(0:NIfTot), ilut_temp(0:NIfTot)
        INTEGER J, JR, ExcitLevel, iunit
        INTEGER LSCR, LISCR, MaxIndex
        LOGICAL tMC!,TestClosedShellDet,Found,tSign
        real(dp) :: calcmcen, calcdlwdb, norm, temp_hel
        integer:: ic, TempnI(NEl), MomSymDet(NEl), ICSym, ICConnect, PairedUnit, SelfInvUnit
        integer(n_int) :: iLutMomSym(0:NIfTot)
        logical :: tSuccess
        type(DavidsonCalcType) :: davidsonCalc
        integer(n_int), allocatable, dimension(:, :) :: davidson_ilut
        integer, allocatable, dimension(:) :: davidson_parities
        integer :: nI(nel)
        integer :: davidson_size
        !character (len=*), parameter :: t_r = "DoDetCalc"
        !integer(TagIntType) :: IlutTag

        IF (tEnergy) THEN
            write(stdout, '(1X,A,E19.3)') ' B2LIMIT : ', B2L
            write(stdout, *) ' NBLK : ', NBLK
            write(stdout, *) ' NKRY : ', NKRY
            write(stdout, *) ' NEVAL : ', NEVAL

            write(stdout, *) ' NCYCLE : ', NCYCLE
            write(stdout, *) ' TENERGY : ', TENERGY
            write(stdout, *) ' IOBS : ', IOBS
            write(stdout, *) ' JOBS : ', JOBS
            write(stdout, *) ' KOBS : ', KOBS
            write(stdout, *) ' NMSH : ', NMSH
            IF (IOBS > NMSH .OR. IOBS <= 0 .OR. JOBS > NMSH .OR. JOBS <= 0 .OR. KOBS > NMSH .OR. KOBS <= 0) THEN
                call stop_all(this_routine, ' !!! REFERENCE PARTICLE NOT IN BOX !!! ')
            end if
        end if

!C.. now back to the storing H
        IF (TCALCHMAT) THEN
            if (t_new_real_space_hubbard) then
                call init_real_space_hubbard()
            end if
            write(stdout, *) "Calculating H matrix"
!C..We need to measure HAMIL and LAB first
            allocate(NROW(NDET), stat=ierr)
            CALL LogMemAlloc('NROW', NDET, 4, this_routine, NROWTag, ierr)
            NROW(1:NDET) = 0
            ICMAX = 1
!Falsify tMC
            TMC = .FALSE.
            allocate(HAMIL(0), LAB(0))
            CALL DETHAM(NDET, NEL, NMRKS, HAMIL, LAB, NROW, .TRUE., ICMAX, GC, TMC)
            deallocate(HAMIL, LAB)
            write(stdout, *) ' FINISHED COUNTING '
            write(stdout, *) "Allocating memory for hamiltonian: ", GC * 2
            CALL neci_flush(stdout)
!C..Now we know size, allocate memory to HAMIL and LAB
            LENHAMIL = GC
            allocate(Hamil(LenHamil), stat=ierr)
            LogAlloc(ierr, 'HAMIL', LenHamil, HElement_t_sizeB, tagHamil)
            HAMIL = (0.0_dp)
!C..
            allocate(LAB(LENHAMIL), stat=ierr)
            CALL LogMemAlloc('LAB', LenHamil, 4, this_routine, LabTag, ierr)

            LAB(1:LENHAMIL) = 0
!C..Now we store HAMIL and LAB
            CALL DETHAM(NDET, NEL, NMRKS, HAMIL, LAB, NROW, .FALSE., ICMAX, GC, TMC)

            IF (BTEST(ILOGGING, 7)) THEN
!C.. we write out H now
                iunit = get_free_unit()
                open(iunit, FILE='HAMIL', STATUS='UNKNOWN')
                J = 0
                JR = 0
!C            HMAX=-dflmax()
!C            HMIN=dflmax()
                DO I = 1, LENHAMIL
                    DO WHILE (I > J)
                        JR = JR + 1
                        J = J + NROW(JR)
                    end do
                    write(iunit, "(2I12)", advance='no') JR, LAB(I)
                    IF (HElement_t_size == 1) THEN
                        write(iunit, *) HAMIL(I)
                    ELSE
                        write(iunit, *) HAMIL(I), ABS(HAMIL(I))
                    end if
!C               CALL WRITEDET(14,NMRKS(1,JR),NEL,.FALSE.)
!C               write(14,"(A)",advance='no'),"|"
!C               CALL WRITEDET(14,NMRKS(1,LAB(I)),NEL,.FALSE.)
!C              write(14,"(F27.20)") HAMIL(I)
!C               CALL WRITEDET(14,NMRKS(1,LAB(I)),NEL,.FALSE.)
!C               write(14,"(A)",advance='no'),"|"
!C               CALL WRITEDET(14,NMRKS(1,JR),NEL,.FALSE.)
!C               write(14,"(F27.20)") HAMIL(I)

!C               IF(HAMIL(I).GT.HMAX) HMAX=HAMIL(I)
!C               IF(HAMIL(I).LT.HMIN) HMIN=HAMIL(I)
                end do
                close(iunit)
            end if
            temp_hel = real(GETHELEMENT(IFDET, IFDET, HAMIL, LAB, NROW, NDET), dp)
            write(stdout, *) '<D0|H|D0>=', temp_hel
            write(stdout, *) '<D0|T|D0>=', CALCT(NMRKS(1 : 1 + nEl, IFDET), NEL)
            CALL neci_flush(stdout)
!CC         CALL HAMHIST(HMIN,HMAX,LENHAMIL,NHISTBOXES)
        end if
!C.. We've now finished calculating H if we were going to.
!C.. IF ENERGY CALC (for which we need to have calced H)
        IF (TENERGY) THEN
            IF (NBLK /= 0) THEN
!C..Things needed for Friesner-Pollard diagonalisation
                IF (TMC) call stop_all(this_routine, 'TMC and TENERGY set - Stopping')
                IF (HElement_t_size /= 1) call stop_all(this_routine, 'Cannot do Lanczos on Complex orbitals.')
                NKRY1 = NKRY + 1
                NBLOCK = MIN(NEVAL, NBLK)
                LSCR = MAX(NDET * NEVAL, 8 * NBLOCK * NKRY)
                LISCR = 6 * NBLOCK * NKRY
!C..
!            write(stdout,'(/,/,8X,64(1H*))')
                write(stdout, '(7X," *",62X,"*")')
                write(stdout, '(7X," *",19X,A,18X,"*")') ' LANCZOS DIAGONALISATION '
                write(stdout, '(7X," *",62X,"*")')
!            write(stdout,'(7X,1X,64(1H*))')
!C..Set up memory for FRSBLKH

                allocate(A(NEVAL, NEVAL), stat=ierr)
                CALL LogMemAlloc('A', NEVAL**2, 8, this_routine, ATag, ierr)
                A = 0.0_dp
!C..
!C,, W is now allocated with CK
!C..
                allocate(V(NDET * NBLOCK * NKRY1), stat=ierr)
                CALL LogMemAlloc('V', NDET * NBLOCK * NKRY1, 8, this_routine, VTag, ierr)
                V = 0.0_dp
!C..
                allocate(AM(NBLOCK * NBLOCK * NKRY1), stat=ierr)
                CALL LogMemAlloc('AM', NBLOCK * NBLOCK * NKRY1, 8, this_routine, AMTag, ierr)
                AM = 0.0_dp
!C..
                allocate(BM(NBLOCK * NBLOCK * NKRY), stat=ierr)
                CALL LogMemAlloc('BM', NBLOCK * NBLOCK * NKRY, 8, this_routine, BMTag, ierr)
                BM = 0.0_dp
!C..
                allocate(T(3 * NBLOCK * NKRY * NBLOCK * NKRY), stat=ierr)
                CALL LogMemAlloc('T', 3 * NBLOCK * NKRY * NBLOCK * NKRY, 8, this_routine, TTag, ierr)
                T = 0.0_dp
!C..
                allocate(WT(NBLOCK * NKRY), stat=ierr)
                CALL LogMemAlloc('WT', NBLOCK * NKRY, 8, this_routine, WTTag, ierr)
                WT = 0.0_dp
!C..
                allocate(SCR(LScr), stat=ierr)
                CALL LogMemAlloc('SCR', LScr, 8, this_routine, SCRTag, ierr)
                SCR = 0.0_dp
                allocate(ISCR(LIScr), stat=ierr)
                CALL LogMemAlloc('IScr', LIScr, 4, this_routine, IScrTag, ierr)
                ISCR(1:LISCR) = 0
                allocate(INDEX(NEVAL), stat=ierr)
                CALL LogMemAlloc('INDEX', NEVAL, 4, this_routine, INDEXTag, ierr)
                INDEX(1:NEVAL) = 0
!C..
                allocate(WH(NDET), stat=ierr)
                CALL LogMemAlloc('WH', NDET, 8, this_routine, WHTag, ierr)
                WH = 0.0_dp
                allocate(WORK2(3 * NDET), stat=ierr)
                CALL LogMemAlloc('WORK2', 3 * NDET, 8, this_routine, WORK2Tag, ierr)
                WORK2 = 0.0_dp
                allocate(V2(NDET, NEVAL), stat=ierr)
                CALL LogMemAlloc('V2', NDET * NEVAL, 8, this_routine, V2Tag, ierr)
                V2 = 0.0_dp
!C..Lanczos iterative diagonalising routine
                if (t_non_hermitian_2_body) then
                    call stop_all(this_routine, &
                                  "NECI_FRSBLKH not adapted for non-hermitian Hamiltonians")
                end if
#ifdef CMPLX_
            call stop_all(this_routine, "this does not make sense for complex code")
#else
                CALL NECI_FRSBLKH(NDET, ICMAX, NEVAL, HAMIL, LAB, CK, CKN, NKRY, NKRY1, NBLOCK, NROW, LSCR, LISCR, A, W, V, AM, BM, T, WT, &
         &  SCR, ISCR, INDEX, NCYCLE, B2L, .true., .false., .true.)
#endif

!Multiply all eigenvalues by -1.
                CALL DSCAL(NEVAL, -1.0_dp, W, 1)
            ELSE
!C.. We splice in a non-Lanczos diagonalisin routine if NBLOCK=0
                IF (NEVAL /= NDET) THEN
                    write(stdout, *) "NEVAL.NE.NDET.", NEVAL, NDET, " Cannot exactly diagonalize."
                    call stop_all(this_routine, "Cannot exactly diagonalise")
                end if
                write(stdout, *) "NBLK=0.  Doing exact diagonalization."
                IF (TCALCHMAT) THEN
                    allocate(WORK(4 * NDET), stat=ierr)
                    CALL LogMemAlloc('WORK', 4 * NDET, 8 * HElement_t_size, this_routine, WorkTag, ierr)
                    allocate(WORK2(3 * NDET), stat=ierr)
                    CALL LogMemAlloc('WORK2', 3 * NDET, 8, this_routine, WORK2Tag, ierr)
                    if (t_non_hermitian_2_body) then
                        call stop_all(this_routine, &
                                      "HDIAG_nec is not setup for non-hermitian Hamiltonians")
                    end if
                    CALL HDIAG_neci(NDET, HAMIL, LAB, NROW, CK, W, WORK2, WORK, NBLOCKSTARTS, NBLOCKS)
                end if
            end if
!C..
!  Since we no longer use HAMIL or LAB, we deallocate
            if (.not. tCalcVariationalEnergy) then
                LogDealloc(tagHamil)
                Deallocate(Hamil)
                DEallocate(LAB)
                CALL LogMemDealloc(this_routine, LabTag)
            end if
            allocate(TKE(NEVAL), stat=ierr)
            CALL LogMemAlloc('TKE', NEVAL, 8, this_routine, TKETag, ierr)

            EXEN = CALCMCEN(NEVAL, W, BETA)
            write(stdout, "(A,F19.9)") "EXACT E(BETA)=", EXEN
            GSEN = CALCDLWDB(IFDET, NDET, NEVAL, CK, W, BETA)
            write(stdout, "(A,F19.9)") "EXACT DLWDB(D0)=", GSEN
            write(stdout, "(A,F19.9)") "GROUND E=", W(1)
!C.. END ENERGY CALC
!      end if
        else if (tFCIDavidson) THEN
            if (.not. TSPN .or. LMS /= 0) then
                call stop_all("DoDetCalc", "FCI-Davidson only works for closed shell systems.")
            end if
            davidsonCalc = davidson_direct_ci_init()
            davidson_size = davidsonCalc%super%space_size
            allocate(davidson_ilut(0:NIfTot, davidson_size))
            allocate(davidson_parities(davidson_size))
            call generate_entire_ras_space(davidson_ras, davidson_classes, davidson_size, davidson_ilut, davidson_parities)
            !Find HF index
            !Set this flag, otherwise hfindex will be overwritten
            tCalcHFIndex = .false.
            davidsonCalc%super%hfindex = 0
            CALL EncodeBitDet(FDet, iLut(0:nifd))
            do i = 1, davidson_size
                if (DetBitEq(davidson_ilut(:, i), ilut)) then
                    davidsonCalc%super%hfindex = i
                    exit
                end if
            end do
            IF (davidsonCalc%super%hfindex == 0) call stop_all("DoDetCalc", "Fermi determinant is not found in RAS space!")
            if (t_non_hermitian_2_body) then
                call stop_all(this_routine, &
                              "perform_davidson not adapted for non-hermitian Hamiltonians!")
            end if
            call perform_davidson(davidsonCalc, direct_ci_type, .true.)
        end if

        call neci_flush(stdout)
!C.. If we're calculating rhos (for which we have to have calced H
!No longer used
!      IF(TRHOIJ) THEN
!         IF((.NOT.TENERGY).AND.(.NOT.TREADRHO)) THEN
!            write(stdout,*) "Calculating approx RHOs"
!            write(stdout,*) "Using Trotter decomposition? ",TTROT
!            write(stdout,*) "Order of Taylor: ",ABS(NTAY)
!            CALL CALCAPPROXRHOIJ(BETA,I_P,HAMIL,LAB,NROW,NDET,RHOMIN,RHOMAX,NRHOS,RHOEPS,TTROT,NTAY)
!         end if
!      end if
!C..Free HAMIL AND LAB memory if we no longer need them
!      IF(TCALCHMAT.AND..NOT.(TMONTE.AND.TMC)) THEN
!      end if

!C.. IF we want to compress the found determinants for use later...
        IF (tFindDets) THEN
            IF (tCompressDets) THEN
!Need to find symmetry of the reference determinant, so that we can only look for determinants of the correct symmetry.
                CALL GETSYM(FDET, NEL, G1, NBASISMAX, IHFSYM)
                IF (.not. associated(NMRKS)) THEN
                    write(stdout, *) "NMRKS not allocated"
                    CALL neci_flush(stdout)
                    CALL Stop_All("DoDetCalc", "NMRKS not allocated so cannot compress dets.")
                end if
!First, we want to count the number of determinants of the correct symmetry...
                Det = 0
                norm = 0.0_dp
                if (tFCIDavidson) then
                    do i = 1, davidson_size
                        CALL decode_bit_det(nI, davidson_ilut(:, i))
                        CALL GETSYM(nI, NEL, G1, NBASISMAX, ISYM)
                        !CALL GetLz(nI,NEL,Lz)
                        IF (ISym%Sym%S == IHFSYM%Sym%S) THEN
                            Det = Det + 1
                            norm = norm + (davidsonCalc%davidson_eigenvector(i))**2
                        end if
                    end do
                else
                    do i = 1, NDET
                        CALL GETSYM(NMRKS(:, i), NEL, G1, NBASISMAX, ISYM)
                        IF (ISym%Sym%S == IHFSYM%Sym%S) THEN
                            Det = Det + 1
                            IF (tEnergy) THEN
                                norm = norm + (REAL(CK(i, 1), dp))**2
                            end if
                        end if
                    end do
                end if
                write(stdout, "(I25,A,I4,A)") Det, " determinants of symmetry ", IHFSym%Sym%S, " found."
                write(stdout, *) "Normalization of eigenvector 1 is: ", norm
                CALL neci_flush(stdout)

                allocate(FCIDets(0:NIfTot, Det), stat=ierr)
                IF (ierr /= 0) CALL Stop_All("DetCalc", "Cannot allocate memory to hold vector")
                allocate(FCIDetIndex(0:NEl + 1), stat=ierr) !+1 so we can store the end of the array too
                allocate(Temp(Det), stat=ierr)
                IF (ierr /= 0) CALL Stop_All("DetCalc", "Cannot allocate memory to hold vector")
                IF (tEnergy .or. tFCIDavidson) THEN
                    allocate(FCIGS(Det), stat=ierr)
                    IF (ierr /= 0) CALL Stop_All("DetCalc", "Cannot allocate memory to hold vector")
                end if
                if (tCalcVariationalEnergy) then
                    !This allows us to resort to get back to the hamiltonian ordering
                    allocate(ReIndex(Det), stat=ierr)
                    if (ierr /= 0) CALL Stop_All("DetCalc", "Cannot allocate memory to hold vector")
                    ReIndex(:) = 0
                end if

                Det = 0
                FCIDetIndex(:) = 0
                if (tFCIDavidson) then
                    do i = 1, davidson_size
                        CALL decode_bit_det(nI, davidson_ilut(:, i))
                        CALL GETSYM(nI, NEL, G1, NBASISMAX, ISYM)
                        !CALL GetLz(nI,NEL,Lz)
                        !IF((ISym%Sym%S.eq.IHFSYM%Sym%S).and.(Lz.eq.LzTot)) THEN
                        IF (ISym%Sym%S == IHFSYM%Sym%S) THEN
                            Det = Det + 1
                            ExcitLevel = FindBitExcitLevel(davidson_ilut(:, i), iLut) !iGetExcitLevel_2(FDet,nI,NEl,NEl)
                            FCIDetIndex(ExcitLevel + 1) = FCIDetIndex(ExcitLevel + 1) + 1
                            Temp(Det) = ExcitLevel    !Temp will now temporarily hold the excitation level of the determinant.
                            FCIDets(:, Det) = davidson_ilut(:, i)
                            FCIGS(Det) = davidson_parities(i) * davidsonCalc%davidson_eigenvector(i) / norm
                        end if
                        if (tCalcVariationalEnergy) ReIndex(i) = i
                    end do
                else
                    do i = 1, NDet
                        CALL GETSYM(NMRKS(:, i), NEL, G1, NBASISMAX, ISYM)
                        IF (ISym%Sym%S == IHFSYM%Sym%S) THEN
                            Det = Det + 1
                            ExcitLevel = iGetExcitLevel_2(FDet, NMRKS(:, i), NEl, NEl)
                            ! FCIDetIndex is off by one, for later cumulative indexing
                            FCIDetIndex(ExcitLevel + 1) = FCIDetIndex(ExcitLevel + 1) + 1
                            Temp(Det) = ExcitLevel    !Temp will now temporarily hold the excitation level of the determinant.
                            CALL EncodeBitDet(NMRKS(:, i), FCIDets(0:NIfTot, Det))
                            IF (tEnergy) THEN
                                FCIGS(Det) = REAL(CK(i, 1), dp) / norm
                            end if
                        end if
                        if (tCalcVariationalEnergy) ReIndex(i) = i
                    end do
                end if

                IF (iExcitLevel <= 0) THEN
                    MaxIndex = NEl
                ELSE
                    MaxIndex = MIN(iExcitLevel, NEl)
                end if
!NB FCIDetIndex is off by 1 for later cumulation
                do i = 1, MaxIndex
                    write(stdout, *) "Number at excitation level: ", i, " is: ", FCIDetIndex(i + 1)
                end do

                ! We now want to sort the determinants according to the
                ! excitation level (stored in Temp)
                IF (.not. tEnergy .and. .not. tFCIDavidson) THEN
                    call sort(temp(1:Det), FCIDets(:, 1:Det))
                ELSE
                    if (tCalcVariationalEnergy) then
                        call sort(temp(1:Det), FCIDets(:, 1:Det), FCIGS(1:Det), ReIndex(1:Det))
                    else
                        call sort(temp(1:Det), FCIDets(:, 1:Det), FCIGS(1:Det))
                    end if
                end if

!Test that HF determinant is the first determinant
                CALL EncodeBitDet(FDet, iLut(0:nifd))
                do i = 0, nifd
                    IF (iLut(i) /= FCIDets(i, 1)) THEN
                        CALL Stop_All("DetCalc", "Problem with ordering the determinants by excitation level")
                    end if
                end do

!Change it so that FCIDetIndex indexes the start of the block of that excitation level.
                FCIDetIndex(0) = 1
                do i = 1, NEl + 1
                    FCIDetIndex(i) = FCIDetIndex(i - 1) + FCIDetIndex(i)
                end do
                FCIDetIndex(nEl + 1) = Det + 1
                IF (FCIDetIndex(nEl + 1) /= Det + 1) THEN
                    write(stdout, *) "FCIDETIndex:", FCIDetIndex(:)
                    CALL Stop_All("DetCalc", "Error in the indexing of determinant excitation level")
                end if

!We now need to sort within the excitation level by the "number" of the determinant
                do i = 1, MaxIndex
                    IF (.not. tEnergy .and. .not. tFCIDavidson) THEN
                        call sort(FCIDets(:, FCIDetIndex(i):FCIDetIndex(i + 1) - 1), &
                                  temp(FCIDetIndex(i):FCIDetIndex(i + 1) - 1))
                    ELSE
                        if (tCalcVariationalEnergy) then
                            call sort(FCIDets(:, FCIDetIndex(i):FCIDetIndex(i + 1) - 1), &
                                      temp(FCIDetIndex(i):FCIDetIndex(i + 1) - 1), &
                                      FCIGS(FCIDetIndex(i):FCIDetIndex(i + 1) - 1), &
                                      ReIndex(FCIDetIndex(i):FCIDetIndex(i + 1) - 1))
                        else
                            call sort(FCIDets(:, FCIDetIndex(i):FCIDetIndex(i + 1) - 1), &
                                      temp(FCIDetIndex(i):FCIDetIndex(i + 1) - 1), &
                                      FCIGS(FCIDetIndex(i):FCIDetIndex(i + 1) - 1))
                        end if
                    end if
                end do

                IF (tEnergy .or. tFCIDavidson) THEN
                    IF (tLogDETS .and. iProcIndex == 0) THEN
                        iunit = get_free_unit()
                        open(iunit, FILE='SymDETS', STATUS='UNKNOWN')

                        do i = 1, Det
                            write(iunit, "(2I17)", advance='no') i, temp(i)
                            do j = 0, nifd
                                write(iunit, "(I17)", advance='no') FCIDets(j, i)
                            end do
                            write(iunit, "(A,G25.16,A)", advance='no') " ", FCIGS(i), "  "
                            Call WriteBitDet(iunit, FCIDets(:, i), .true.)
                        end do
                        close(iunit)
                    end if
                    DEallocate(FCIGS)
                ELSE
                    IF (tLogDETS .and. iProcIndex == 0) THEN
                        iunit = get_free_unit()
                        open(iunit, FILE='SymDETS', STATUS='UNKNOWN')
                        write(iunit, *) "FCIDETIndex: ", FCIDetIndex(:)
                        write(iunit, *) "***"
                        do i = 1, Det
                            write(iunit, "(2I13)", advance='no') i, temp(i)
                            do j = 0, nifd
                                write(iunit, "(I13)", advance='no') FCIDets(j, i)
                            end do
                            write(iunit, "(A)", advance='no') " "
                            Call WriteBitDet(iunit, FCIDets(:, i), .true.)
                        end do
                        close(iunit)
                    end if
                end if !tEnergy - for dumping compressed ordered GS wavefunction
                DEallocate(Temp)
            end if !tCompressDets
        end if !tFindDets
!C..
        IF (TEnergy) THEN
            CALL CFF_CHCK(NDET, NEVAL, NMRKS, NEL, G1, CK, TKE)
            IF (BTEST(ILOGGING, 7)) CALL WRITE_PSI(BOX, BOA, COA, NDET, NEVAL, NBASISMAX, NEL, CK, W)
            IF (BTEST(ILOGGING, 8)) CALL WRITE_PSI_COMP(BOX, BOA, COA, NDET, NEVAL, NBASISMAX, NEL, CK, W)
            write(stdout, *) '       ====================================================== '
            write(stdout, '(A5,5X,A15,1X,A18,1x,A20)') 'STATE', 'KINETIC ENERGY', 'COULOMB ENERGY', 'TOTAL ENERGY'
            iunit = get_free_unit()
            open(iunit, FILE='ENERGIES', STATUS='UNKNOWN')
            DO IN = 1, NEVAL
                write(stdout, '(I5,2X,3(F19.11,2x))') IN, TKE(IN), W(IN) - TKE(IN), W(IN)
                !            write(iunit,"(I7)",advance='no') IN
                !            CALL WRITEDET(iunit,NMRKS(1,IN),NEL,.FALSE.)
                write(iunit, "(F19.11)") W(IN)
            end do
            close(iunit)
            write(stdout, *) '       ====================================================== '
        else if (tFCIDavidson) THEN
            open(iunit, FILE='ENERGIES', STATUS='UNKNOWN')
            write(iunit, "(F19.11)") davidsonCalc%davidson_eigenvalue
            close(iunit)
!C., END energy calc
        end if

        IF (tFCIDavidson) THEN
            call davidson_direct_ci_end(davidsonCalc)
            call DestroyDavidsonCalc(davidsonCalc)
            DEallocate(davidson_ilut)
            DEallocate(davidson_parities)
        end if
!C.. Jump to here if just read Psi in
        CONTINUE

        !Now deallocate NMRKS if tFindDets and not tEnergy
        if (tFindDets .and. tCompressDets .and. (.not. tEnergy)) then
            DEallocate(NMRKS)
            CALL LogMemDealloc(this_routine, tagNMRKS)
        end if

    End Subroutine DoDetCalc