FillCoeffT1 Subroutine

    subroutine FillCoeffT1

        use RotateOrbsData, only: CoeffT1, SymLabelList3_rot, SymOrbs_rot, SymOrbs_rotTag, &
                                  TruncEval, NoRotOrbs, EvaluesTrunc, EvaluesTruncTag
        use LoggingData, only: tTruncRODump, tTruncDumpbyVal

        integer :: l, k, i, j, NoRotAlphBet, io1, io2, ierr
        character(len=*), parameter :: t_r = 'FillCoeffT1'
        character(len=5) :: Label
        character(len=20) :: LabelFull
        real(dp) :: OccEnergies(1:NoRotOrbs)

        FillCoeff_Time%timer_name = 'FillCoeff'
        call set_timer(FillCoeff_Time, 30)

        write(stdout, *) 'NatOrbMat'
        do i = 1, nBasis
            do j = 1, nBasis
                write(stdout, '(F10.6)', advance='no') NatOrbMat(j, i)
            end do
            write(stdout, *) ''
        end do

        if (tTruncRODump) then

            if (tTruncDumpbyVal) then
                NoFrozenVirt = 0
                if (tStoreSpinOrbs) then
                    do i = SpatOrbs, 1, -1
                        if (Evalues(i) > TruncEval) exit
                        if (Evalues(i + SpatOrbs) > TruncEval) exit
                        NoFrozenVirt = NoFrozenVirt + 2
                    end do
                    if (NoFrozenVirt >= (NoOrbs - NEl)) call stop_all(t_r, 'Freezing all virtual orbitals.')
                else
                    do i = SpatOrbs, 1, -1
                        if (Evalues(i) > TruncEval) exit
                        NoFrozenVirt = NoFrozenVirt + 1
                    end do
                    if (NoFrozenVirt >= (SpatOrbs - (NEl / 2))) call stop_all(t_r, 'Freezing all virtual orbitals.')
                end if
                NoRotOrbs = NoOrbs - NoFrozenVirt
            end if

            allocate(SymOrbs_rot(NoOrbs), stat=ierr)
            call LogMemAlloc('SymOrbs_rot', NoOrbs, 4, t_r, SymOrbs_rotTag, ierr)
            SymOrbs_rot(:) = 0

            allocate(EvaluesTrunc(NoOrbs - NoFrozenVirt), stat=ierr)
            call LogMemAlloc('EvaluesTrunc', NoOrbs - NoFrozenVirt, 4, t_r, EvaluesTruncTag, ierr)
            EvaluesTrunc(:) = 0.0_dp

            if (tStoreSpinOrbs) then
                NoRotAlphBet = SpatOrbs - (NoFrozenVirt / 2)
            else
                NoRotAlphBet = NoOrbs - NoFrozenVirt
            end if

            if (tStoreSpinOrbs) then
                ! As we reorder these so that they are truncated, we also need
                ! to pair up symmetries.

                call CalcOccEnergies(OccEnergies)

                ! First nOccBeta, then nOccAlpha.
                do i = 1, (2 * nOccBeta), 2
                    k = 1
                    do while (abs(OccEnergies(k)) < 1.0e-10_dp)
                        k = k + 2
                    end do
                    do j = 1, (2 * nOccBeta), 2
                        if ((OccEnergies(j) < OccEnergies(k)) .and. (abs(OccEnergies(j)) > 1.0e-10_dp)) k = j
                    end do
                    l = ceiling(real(k, dp) / 2.0_dp)
                    CoeffT1(:, i) = NatOrbMat(:, l)
                    EvaluesTrunc(i) = Evalues(l)
                    SymOrbs_rot(i) = SymOrbs_rotTemp(l)
                    OccEnergies(k) = 0.0_dp
                end do
                do i = 2, (2 * nOccAlpha), 2
                    k = 2
                    do while (abs(OccEnergies(k)) < 1.0e-10_dp)
                        k = k + 2
                    end do
                    do j = 2, (2 * nOccAlpha), 2
                        if ((OccEnergies(j) < OccEnergies(k)) .and. (abs(OccEnergies(j)) > 1.0e-10_dp)) k = j
                    end do
                    l = (k / 2) + SpatOrbs
                    CoeffT1(:, i) = NatOrbMat(:, l)
                    EvaluesTrunc(i) = Evalues(l)
                    SymOrbs_rot(i) = SymOrbs_rotTemp(l)
                    OccEnergies(k) = 0.0_dp
                end do

                ! Need to fill coeffT1 so that it goes alpha beta alpha beta.
                k = (2 * nOccBeta) + 1
                do i = nOccBeta + 1, NoRotAlphBet
                    CoeffT1(:, k) = NatOrbMat(:, i)
                    EvaluesTrunc(k) = Evalues(i)
                    SymOrbs_rot(k) = SymOrbs_rotTemp(i)
                    k = k + 2
                end do
                k = (2 * nOccAlpha) + 2
                do i = SpatOrbs + nOccAlpha + 1, SpatOrbs + NoRotAlphBet
                    CoeffT1(:, k) = NatOrbMat(:, i)
                    SymOrbs_rot(k) = SymOrbs_rotTemp(i)
                    EvaluesTrunc(k) = Evalues(i)
                    k = k + 2
                end do

            else

                ! Order occupied in terms of energy again - this makes sure
                ! freezing etc doesn't get screwed up.
                call CalcOccEnergies(OccEnergies)

                ! OccEnergies has the orbital energies as they are ordered
                ! currently - need to put NatOrbMat into CoeffT1 so that this
                ! goes from lowest energy to highest.

                do i = 1, NEl / 2
                    k = 1
                    do while (abs(OccEnergies(k)) < 1.0e-10_dp)
                        k = k + 1
                    end do
                    do j = 1, NEl / 2
                        if ((OccEnergies(j) < OccEnergies(k)) .and. (abs(OccEnergies(j)) > 1.0e-10_dp)) k = j
                    end do
                    CoeffT1(:, i) = NatOrbMat(:, k)
                    EvaluesTrunc(i) = Evalues(k)
                    SymOrbs_rot(i) = SymOrbs_rotTemp(k)
                    OccEnergies(k) = 0.0_dp
                end do

                do i = (NEl / 2) + 1, NoRotAlphBet
                    CoeffT1(:, i) = NatOrbMat(:, i)
                    EvaluesTrunc(i) = Evalues(i)
                    SymOrbs_rot(i) = SymOrbs_rotTemp(i)
                end do
            end if

        else

            do i = 1, NoOrbs
                CoeffT1(:, i) = NatOrbMat(:, i)
            end do

        end if

        if (tTruncRODump) then

            Label = ''
            LabelFull = ''
            write(Label, '(I5)') NoFrozenVirt
            LabelFull = 'EVALUES-TRUNC-'//adjustl(Label)

            io1 = get_free_unit()
            open(io1, file=LabelFull, status='unknown')
            if (tStoreSpinOrbs) then
                write(io1, *) NoOrbs - NoFrozenVirt
                do i = 1, NoOrbs - NoFrozenVirt, 2
                    write(io1, '(I5,ES20.10,I5,A5,I5,ES20.10,I5)') i, EvaluesTrunc(i), SymOrbs_rot(i), &
                        '  *  ', i + 1, EvaluesTrunc(i + 1), SymOrbs_rot(i + 1)
                end do
            else
                write(io1, *) NoOrbs - NoFrozenVirt
                do i = 1, NoOrbs - NoFrozenVirt
                    write(io1, '(ES20.10,I5)') EvaluesTrunc(i), SymOrbs_rot(i)
                end do
            end if
            close(io1)
        else
            io2 = get_free_unit()
            open(io2, file='EVALUES', status='unknown')
            write(io2, *) NoOrbs
            if (tStoreSpinOrbs) then
                k = 0
                do i = 1, NoOrbs, 2
                    k = k + 1
                    if (tTruncRODump) then
                        write(io2, '(2I5,ES20.10,I5,A5,I5,ES20.10,I5)') (NoOrbs - i + 1), i, Evalues(k), SymOrbs_rot(i), &
                            '  *  ', i + 1, Evalues(k + SpatOrbs), SymOrbs_rot(i + 1)
                    else
                        write(io2, '(2I5,ES20.10,I5,A5,I5,ES20.10,I5)') (NoOrbs - i + 1), i, Evalues(k), &
                            int(G1(SymLabelList3_rot(k))%Sym%S, 4), '  *  ', &
                            i + 1, Evalues(k + SpatOrbs), int(G1(SymLabelList3_rot(k + SpatOrbs))%Sym%S, 4)
                    end if
                end do
            else
                do i = 1, SpatOrbs
                    write(io2, '(3I5,ES20.10)') i, NoOrbs - i + 1, (NoOrbs - i + 1) * 2, Evalues(i)
                end do
            end if
            close(io2)
        end if

        call HistNatOrbEvalues()

        call halt_timer(FillCoeff_Time)

    end subroutine FillCoeffT1