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 (non_hermitian()) 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 (non_hermitian()) 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 (non_hermitian()) 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