DiagNatOrbMat Subroutine

    subroutine DiagNatOrbMat()

        ! The diagonalisation routine reorders the orbitals in such a way that
        ! the corresponding orbital labels are lost. In order to keep the spin
        ! and spatial symmetries, each symmetry must be fed into the
        ! diagonalisation routine separately. The best way to do this is to
        ! order the orbitals so that all the alpha orbitals follow all the beta
        ! orbitals, with the occupied orbitals first, in terms of symmetry, and
        ! the virtual second, also ordered by symmetry. This gives us
        ! flexibility w.r.t rotating only the occupied or only virtual and
        ! looking at high spin states.

        use MemoryManager, only: TagIntType

        real(dp) :: SumTrace, SumDiagTrace
        real(dp), allocatable :: WORK2(:), EvaluesSym(:), NOMSym(:, :)
        integer :: ierr, i, j, x, z, Sym, LWORK2
        integer :: SymStartInd, NoSymBlock, PrevSym, StartOccVirt, EndOccVirt, Prev, NoOcc
        integer(TagIntType) :: EvaluesSymTag, NOMSymTag, WORK2Tag
        character(len=*), parameter :: t_r = 'DiagNatOrbMat'

        DiagNatOrbMat_Time%timer_name = 'DiagNatOrb'
        call set_timer(DiagNatOrbMat_Time, 30)

        do x = 1, NoSpinCyc
            if (tSeparateOccVirt) then
                if (x == 1) then
                    if (tStoreSpinOrbs) then
                        NoOcc = nOccBeta
                    else
                        NoOcc = NEl / 2
                    end if
                    Prev = 0
                else if (x == 2) then
                    NoOcc = nOccAlpha
                    Prev = SpatOrbs
                end if
            else
                NoOcc = 0
            end if
            if (tRotateVirtOnly) then
                do i = 1, NoOcc
                    do j = 1, SpatOrbs
                        NatOrbMat(i + Prev, j + Prev) = 0.0_dp
                        NatOrbMat(j + Prev, i + Prev) = 0.0_dp
                        if (i == j) NatOrbMat(i + Prev, j + Prev) = 1.0_dp
                    end do
                    Evalues(i + Prev) = 1.0_dp
                end do
            else if (tRotateOccOnly) then
                do i = NoOcc + 1, SpatOrbs
                    do j = 1, SpatOrbs
                        NatOrbMat(i + Prev, j + Prev) = 0.0_dp
                        NatOrbMat(j + Prev, i + Prev) = 0.0_dp
                        if (i == j) NatOrbMat(i + Prev, j + Prev) = 1.0_dp
                    end do
                    Evalues(i + Prev) = 1.0_dp
                end do
            else if (tSeparateOccVirt) then
                do i = 1, NoOcc
                    do j = NoOcc + 1, SpatOrbs
                        NatOrbMat(i + Prev, j + Prev) = 0.0_dp
                        NatOrbMat(j + Prev, i + Prev) = 0.0_dp
                    end do
                end do
            end if
        end do

        ! Test that we're not breaking symmetry.
        do i = 1, NoOrbs
            do j = 1, NoOrbs
                if (tStoreSpinOrbs) then
                    if ((int(G1(SymLabelList2_rot(i))%sym%S, 4) /= int(G1(SymLabelList2_rot(j))%sym%S, 4))) then
                        if (abs(NatOrbMat(i, j)) >= 1.0e-15_dp) then
                            write(stdout, '(6A8,A20)') 'i', 'j', 'Label i', 'Label j', 'Sym i', 'Sym j', 'Matrix value'
                            write(stdout, '(6I3,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
                                int(G1(SymLabelList2_rot(i))%sym%S, 4), &
                                int(G1(SymLabelList2_rot(j))%sym%S, 4), NatOrbMat(i, j)
                            if (tUseMP2VarDenMat) then
                                write(stdout, *) '**WARNING** - There is a non-zero NatOrbMat value between " &
                                 & //"orbitals of different symmetry.'
                                write(stdout, *) 'These elements will be ignored, and the symmetry maintained " &
                                 & //"in the final transformation matrix.'
                            else
                                call stop_all(t_r, 'Non-zero NatOrbMat value between different symmetries.')
                            end if
                        end if
                        NatOrbMat(i, j) = 0.0_dp
                    end if
                else
                    if ((int(G1(SymLabelList2_rot(i) * 2)%sym%S, 4) /= int(G1(SymLabelList2_rot(j) * 2)%sym%S, 4))) then
                        if (abs(NatOrbMat(i, j)) >= 1.0e-15_dp) then
                            write(stdout, '(6A8,A20)') 'i', 'j', 'Label i', 'Label j', 'Sym i', 'Sym j', 'Matrix value'
                            write(stdout, '(6I3,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
                                int(G1(SymLabelList2_rot(i) * 2)%sym%S, 4), int(G1(SymLabelList2_rot(j) * 2)%sym%S, 4), NatOrbMat(i, j)
                            if (tUseMP2VarDenMat) then
                                write(stdout, *) '**WARNING** - There is a non-zero NatOrbMat value between orbitals of " &
                                & //"different symmetry.'
                                write(stdout, *) 'These elements will be ignored, and the symmetry maintained in the " &
                                & //"final transformation matrix.'
                            else
                                call stop_all(t_r, 'Non-zero NatOrbMat value between different symmetries.')
                            end if
                        end if
                        NatOrbMat(i, j) = 0.0_dp
                    end if
                end if
            end do
        end do

        SumTrace = 0.0_dp
        do i = 1, NoOrbs
            SumTrace = SumTrace + NatOrbMat(i, i)
        end do

        write(stdout, *) 'Calculating eigenvectors and eigenvalues of NatOrbMat'
        call neci_flush(stdout)

        ! If we are using spin orbitals, need to feed in the alpha and beta
        ! spins separately. Otherwise these jumble up and the final ordering
        ! is uncorrect. There should be no non-zero values between these, but
        ! can put a check in for this.

        do x = 1, NoSpinCyc

            ! If we want to maintain the symmetry, we cannot have all the
            ! orbitals jumbled up when the diagonaliser reorders the
            ! eigenvectors. Must instead feed each symmetry block in
            ! separately. This means that although the transformed orbitals
            ! are jumbled within the symmetry blocks, the symmetry labels are
            ! all that are relevant and these are unaffected.
            StartOccVirt = 1
            EndOccVirt = 2
            if (tRotateVirtOnly) StartOccVirt = 2
            if (tRotateOccOnly) EndOccVirt = 1

            do z = StartOccVirt, EndOccVirt
                if (x == 1) then
                    if (z == 1) PrevSym = 1
                    if (z == 2) PrevSym = 9
                else if (x == 2) then
                    if (z == 1) PrevSym = 17
                    if (z == 2) PrevSym = 25
                end if

                Sym = 0
                LWORK2 = -1

                do while (Sym <= 7)

                    NoSymBlock = SymLabelCounts2_rot(2, Sym + PrevSym)

                    ! This is one less than the index that the symmetry starts,
                    ! so that when we run through i = 1,...,  we can start at
                    ! SymStartInd+i.
                    SymStartInd = SymLabelCounts2_rot(1, Sym + PrevSym) - 1

                    if (NoSymBlock > 1) then

                        allocate(NOMSym(NoSymBlock, NoSymBlock), stat=ierr)
                        call LogMemAlloc('NOMSym', NoSymBlock**2, 8, t_r, NOMSymTag, ierr)
                        if (ierr /= 0) call stop_all(t_r, "Problem allocating NOMSym.")
                        allocate(EvaluesSym(NoSymBlock), stat=ierr)
                        call LogMemAlloc('EvaluesSym', NoSymBlock, 8, t_r, EvaluesSymTag, ierr)
                        if (ierr /= 0) call stop_all(t_r, "Problem allocating EvaluesSym.")

                        LWORK2 = 3 * NoSymBlock + 1
                        allocate(WORK2(LWORK2), stat=ierr)
                        call LogMemAlloc('WORK2', LWORK2, 8, t_r, WORK2Tag, ierr)
                        if (ierr /= 0) call stop_all(t_r, "Problem allocating WORK2.")

                        do j = 1, NoSymBlock
                            do i = 1, NoSymBlock
                                NOMSym(i, j) = NatOrbMat(SymStartInd + i, SymStartInd + j)
                            end do
                        end do

                        write(stdout, *) '*****'
                        write(stdout, *) 'Symmetry ', Sym, 'with x ', x, ' and z ', z, ' has ', NoSymBlock, ' orbitals.'
                        write(stdout, *) 'The NatOrbMat for this symmetry block is '
                        do i = 1, NoSymBlock
                            do j = 1, NoSymBlock
                                write(stdout, '(F20.10)', advance='no') NOMSym(j, i)
                            end do
                            write(stdout, *) ''
                        end do

                        ! NOMSym goes in as the original NOMSym, comes out as
                        ! the eigenvectors (Coefficients). EvaluesSym comes out
                        ! as the eigenvalues in ascending order.
                        call dsyev('V', 'L', NoSymBlock, NOMSym, NoSymBlock, EvaluesSym, WORK2, LWORK2, ierr)

                        write(stdout, *) 'After diagonalization, the e-vectors (diagonal elements) of this matrix are,'
                        do i = 1, NoSymBlock
                            write(stdout, '(F20.10)', advance='no') EvaluesSym(i)
                        end do
                        write(stdout, *) ''
                        write(stdout, *) 'These go from orbital,', SymStartInd + 1, ' to ', SymStartInd + NoSymBlock

                        do i = 1, NoSymBlock
                            Evalues(SymStartInd + i) = EvaluesSym(i)
                        end do

                        ! CAREFUL if eigenvalues are put in ascending order,
                        ! this may not be correct, with the labelling system.
                        ! Maybe better to just take coefficients and transform
                        ! TMAT2DRot in transform2elints. A check that comes out
                        ! as diagonal is a check of this routine anyway.

                        write(stdout, *) 'The eigenvectors (coefficients) for symmetry block ', Sym
                        do i = 1, NoSymBlock
                            do j = 1, NoSymBlock
                                write(stdout, '(F20.10)', advance='no') NOMSym(j, i)
                            end do
                            write(stdout, *) ''
                        end do

                        ! Directly fill the coefficient matrix with the
                        ! eigenvectors from the diagonalization.
                        do j = 1, NoSymBlock
                            do i = 1, NoSymBlock
                                NatOrbMat(SymStartInd + i, SymStartInd + j) = NOMSym(i, j)
                            end do
                        end do

                        deallocate(WORK2)
                        call LogMemDealloc(t_r, WORK2Tag)

                        deallocate(NOMSym)
                        call LogMemDealloc(t_r, NOMSymTag)

                        deallocate(EvaluesSym)
                        call LogMemDealloc(t_r, EvaluesSymTag)

                    else if (NoSymBlock == 1) then
                        ! The eigenvalue is the lone value, while the
                        ! eigenvector is 1.

                        Evalues(SymStartInd + 1) = NatOrbMat(SymStartInd + 1, SymStartInd + 1)
                        NatOrbMat(SymStartInd + 1, SymStartInd + 1) = 1.0_dp
                        write(stdout, *) '*****'
                        write(stdout, *) 'Symmetry ', Sym, ' has only one orbital.'
                        write(stdout, *) 'Copying diagonal element,', SymStartInd + 1, 'to NatOrbMat'
                    end if

                    Sym = Sym + 1
                end do
            end do
        end do

        write(stdout, *) 'Matrix diagonalised'
        call neci_flush(stdout)

        SumDiagTrace = 0.0_dp
        do i = 1, NoOrbs
            SumDiagTrace = SumDiagTrace + Evalues(i)
        end do
        if ((abs(SumDiagTrace - SumTrace)) > 10.0_dp) then
            write(stdout, *) 'Sum of diagonal NatOrbMat elements : ', SumTrace
            write(stdout, *) 'Sum of eigenvalues : ', SumDiagTrace
            write(stdout, *) 'WARNING, The trace of the 1RDM matrix before diagonalisation is not equal to that after.'
        end if

        call halt_timer(DiagNatOrbMat_Time)

    end subroutine DiagNatOrbMat