IntFreeze Subroutine

    Subroutine IntFreeze
        use SystemData, only: Brr, CoulDampOrb, ECore, fCoulDampMu
        use SystemData, only: G1, iSpinSkip
        use SystemData, only: nBasis, nEl, arr, nbasismax
        use UMatCache, only: GetUMatSize
        use SymData, only: TwoCycleSymGens
        use MemoryManager, only: TagIntType
        use global_utilities

        character(25), parameter ::this_routine = 'IntFreeze'
!//Locals
        HElement_t(dp), pointer :: UMAT2(:)
        INTEGER(TagIntType) tagUMat2
        INTEGER nOcc
        integer(MPIArg) :: umat2_win
        integer(int64) :: UMATInt
        integer nHG

        nHG = nBasis

        if(NEL + 1 < nBasis) CHEMPOT = (ARR(NEL, 1) + ARR(NEL + 1, 1)) / 2.0_dp
!      write(stdout,*) "Chemical Potential: ",CHEMPOT
        IF(NTFROZEN < 0) THEN
            write(stdout, *) "NTFROZEN<0.  Leaving ", -NTFROZEN, " unfrozen virtuals."
            NTFROZEN = NTFROZEN + nBasis - NEL
        end if
        IF((NFROZEN + NFROZENIN) > NEL) CALL Stop_All("IntFreeze", "Overlap between low energy frozen orbitals &
                                  & and inner frozen occupied orbitals - to many frozen occupied orbitals &
                                  & for the number of electrons.")
        IF((NTFROZEN + NTFROZENIN) > (NBASIS - NEL)) CALL Stop_All("IntFreeze", "Overlap between high energy frozen orbitals &
                                  & and inner frozen virtual orbitals - to many frozen virtual orbitals &
                                  & for the number of unnoccupied orbitals.")
        IF(((NFROZENIN > 0) .or. (NTFROZENIN > 0)) .and. (.not. TwoCycleSymGens)) THEN
            CALL Stop_All("IntFreeze", "TwoCycleSymGens is not true. &
            & The code is only set up to deal with freezing from the inside for molecular &
            & systems with only 8 symmetry irreps.")
        end if
        IF(NFROZEN > 0 .OR. NTFROZEN > 0 .OR. NFROZENIN > 0 .OR. NTFROZENIN > 0) THEN
            write(stdout, '(A)') '======== FREEZING ORBITALS =========='
!!C.. At this point, we transform the UMAT and TMAT into a new UMAT and
!!C.. TMAT and Ecore with the frozen orbitals factored in
!!C..
!!C.. a,b are frozen spinorbitals
!!C.. E'core = Ecore+sum_a t_aa + sum_(a<b) (<ab|ab>-<ab|ba>)
!!C.. t'_ii = t_ii+ sum_a ( <ai|ai> - <ai|ia> )
!!C.. NHG contains the old number of orbitals
!!C.. NBASIS contains the new
            NBASIS = NBASIS - NFROZEN - NTFROZEN - NFROZENIN - NTFROZENIN
!!C.. We need to transform some integrals
            !CALL N_MEMORY(IP_TMAT2,HElement_t_size*(NBASIS)**2,'TMAT2')
            !TMAT2=(0.0_dp)
            IF(NBASISMAX(1, 3) >= 0 .AND. ISPINSKIP /= 0) THEN
                CALL GetUMatSize(nBasis, UMATINT)
                call shared_allocate_mpi(umat2_win, umat2,(/UMatInt/))
                !allocate(UMat2(UMatInt), stat=ierr)
                LogAlloc(ierr, 'UMat2', int(UMatInt), HElement_t_SizeB, tagUMat2)
                UMAT2 = 0.0_dp
            ELSE
!!C.. we don't precompute 4-e integrals, so don't allocate a large UMAT
                call shared_allocate_mpi(umat2_win, umat2,(/1_int64/))
                !allocate(UMat2(1), stat=ierr)
                LogAlloc(ierr, 'UMat2', 1, HElement_t_SizeB, tagUMat2)
            end if
!         CALL N_MEMORY_CHECK()

            write(stdout, *) "Freezing ", NFROZEN, " core orbitals."
            write(stdout, *) "Freezing ", NTFROZEN, " virtual orbitals."
            IF(NFROZENIN /= 0) write(stdout, *) "Freezing ", NFROZENIN, " of the highest energy occupied (inner) orbitals."
            IF(NTFROZENIN /= 0) write(stdout, *) "Freezing ", NTFROZENIN, " of the lowest energy virtual (inner) orbitals."

!At the end of IntFREEZEBASIS, NHG is reset to nBasis - the final number of active orbitals.
          CALL IntFREEZEBASIS(NHG, NBASIS, UMAT, UMAT2, ECORE, G1, NBASISMAX, ISPINSKIP, BRR, NFROZEN, NTFROZEN, NFROZENIN, NTFROZENIN, NEL)
            CALL neci_flush(stdout)
            write(stdout, *) "ECORE now", ECORE
            write(stdout, *) "Number of orbitals remaining: ", NBASIS
            nel_pre_freezing = nel
            NEL = NEL - NFROZEN - NFROZENIN
            NOCC = NEL / 2
!!C.. NEL now only includes active electrons
            write(stdout, *) "Number of active electrons:", NEL

            !CALL N_FREEM(IP_TMAT)
            !IP_TMAT=IP_TMAT2
            !IP_TMAT2=NULL
!!C.. Now we can remove the old UMATRIX, and set the pointer UMAT to point
!!C.. to UMAT2
            LogDealloc(tagUMat)
            call shared_deallocate_mpi(int(umat_win, MPIArg), umat)
            !Deallocate(UMat)
            umat_win = umat2_win
            UMat => UMat2
            nullify(UMat2)
            tagUMat = tagUMat2
            tagUMat2 = 0
            call setup_UMatInd()
            CALL writebasis(stdout, G1, NHG, ARR, BRR)
        end if

        ! Setup the umatel pointers as well
        call init_getumatel_fn_pointers()
        call init_bit_rep()

        IF(COULDAMPORB > 0) THEN
            FCOULDAMPMU = (ARR(COULDAMPORB, 1) + ARR(COULDAMPORB + 1, 1)) / 2
            write(stdout, *) "Setting Coulomb damping mu between orbitals ", ARR(COULDAMPORB, 1), " and ", ARR(COULDAMPORB + 1, 1)
            write(stdout, *) "MU=", FCOULDAMPMU
        end if

    End Subroutine IntFreeze