UHFSCF Subroutine

    SUBROUTINE UHFSCF(NBASIS, G1, BRR, ECORE, HFE, HFBASIS, NHFIT, NEL, MS, &
                      FMAT, DMAT, ODMAT, WORK, NSPINS, NSBASIS, HFES, OFMAT, HFMIX, &
                      EDELTA, CDELTA, TRHF, R1, R2, IHFMETHOD, TREADHF, &
                      FRAND, HFDET, ILOGGING)
        INTEGER NSPINS, NSBASIS
        INTEGER NBASIS
        TYPE(BasisFn) G1(*)
        real(dp) ECORE
        INTEGER BRR(NBASIS)
        HElement_t(dp) HFBASIS(NBASIS, NBASIS)
        real(dp) HFES(NSBASIS, NSPINS), HFE(NBASIS)
        HElement_t(dp) FMAT(NSBASIS, NSBASIS, NSPINS)
        HElement_t(dp) OFMAT(NSBASIS, NSBASIS, NSPINS)
        HElement_t(dp) DMAT(NSBASIS, NSBASIS, NSPINS)
        HElement_t(dp) ODMAT(NSBASIS, NSBASIS, NSPINS), WORK(NBASIS * 3)
        real(dp) HFMIX
        HElement_t(dp) R1(NSBASIS, NSBASIS), R2(NSBASIS, NSBASIS)
        INTEGER NHFIT, NEL
        INTEGER I, J, K, L, ISPN, NELS(NSPINS)
        real(dp) ELAST, EDELTA, ECUR, CDELTA
        INTEGER IHFIT
        INTEGER MS, IHFMETHOD, IRHFB
        real(dp) F
        real(dp) RMSD, FRAND
        LOGICAL BR
        LOGICAL TRHF, TREADHF
        INTEGER HFDET(*)
        INTEGER ILOGGING
        F = HFMIX
!         EDELTA=1.0e-8_dp
        ELAST = 1.D20
        WRITE(stdout, *) "Performing Hartree-Fock SCF diagonalisation..."
        IF (IHFMETHOD == -1) THEN
            WRITE(stdout, *) "Method -1:Rotational Mixing "
        ELSEIF (IHFMETHOD == 0) THEN
            WRITE(stdout, *) "Method 0:Linear Mixing"
        END IF
        WRITE(stdout, *) "HF Mixing", HFMIX
        WRITE(stdout, *) "E Thresh:", EDELTA
        WRITE(stdout, *) "RMSD Thresh:", CDELTA
        IF (NSPINS == 2) THEN
            NELS(2) = (MS + NEL) / 2
            NELS(1) = NEL - NELS(2)
            WRITE(stdout, *) " Beta, Alpha: ", NELS(1), NELS(2)
        ELSE
            NELS(1) = NEL
        END IF
        IF (TREADHF) THEN
            CALL READHFFMAT(NBASIS, FMAT, HFES, G1, NSPINS, NSBASIS, .FALSE.)
        ELSE
            CALL GENHFGUESS(FMAT, NSPINS, NSBASIS, BRR, G1, .FALSE., MS, FRAND, NELS, HFDET)
        END IF
        WRITE(stdout, *) "Iteration   Energy     MSD"
        BR = .TRUE.
        IHFIT = 1
        IRHFB = 0
        IF (TRHF .AND. NSPINS > 1) THEN
        IF (NELS(2) > NELS(1)) THEN
            IRHFB = 2
        ELSE
            IRHFB = 1
        END IF
        END IF
        DO WHILE (BR)
        IF (IRHFB > 0) THEN
            CALL DCOPY(NSBASIS * NSBASIS * HElement_t_size, FMAT(1, 1, IRHFB), 1, FMAT(1, 1, 3 - IRHFB), 1)
            CALL DCOPY(NSBASIS * NSBASIS * HElement_t_size, DMAT(1, 1, IRHFB), 1, DMAT(1, 1, 3 - IRHFB), 1)
        END IF
        CALL DCOPY(NSBASIS * NSBASIS * NSPINS * HElement_t_size, DMAT, 1, ODMAT, 1)
        CALL DCOPY(NSBASIS * NSBASIS * NSPINS * HElement_t_size, FMAT, 1, OFMAT, 1)
!.. Construct the Density Matrix
        CALL GENDMAT(NSPINS, NSBASIS, NELS, FMAT, DMAT, .FALSE.)
!.. See how much our density matrix has changed from last time
        RMSD = 0.0_dp
        DO ISPN = 1, NSPINS
        DO I = 1, NSBASIS
        DO J = 1, NSBASIS
            IF (.NOT. TRHF .OR. TRHF .AND. ISPN == IRHFB) RMSD = RMSD + abs(DMAT(I, J, ISPN) - ODMAT(I, J, ISPN))**2
        END DO
        END DO
        END DO
        RMSD = SQRT(RMSD / (NSBASIS * NSBASIS * NSPINS))
!.. replace our HF orbitals in FMAT with the Fock matrix
!.. FMAT just stores the results and the HF orbitals (currently) in it
!.. are not used in calculating the F matrix
        CALL GENFMAT(FMAT, DMAT, NSBASIS, NSPINS)
        CALL DIAGFMAT(NSPINS, NSBASIS, NELS, FMAT, DMAT, HFES, WORK, ECORE, ECUR)
        IF (IRHFB > 0) THEN
            CALL DCOPY(NSBASIS * NSBASIS, FMAT(1, 1, IRHFB), 1, FMAT(1, 1, 3 - IRHFB), 1)
            CALL DCOPY(NSBASIS, HFES(1, IRHFB), 1, HFES(1, 3 - IRHFB), 1)
        END IF
!.. FMAT now contains HF orbitals again, and ECUR the Fock Energy
!.. eigenvector N is in FMAT(i,N,ISPN), where i is the component of the
!.. vector
        WRITE(stdout, *) IHFIT, ECUR, RMSD
!.. Now add back in some of our original F matrix
        IF (IHFMETHOD == -1) THEN
            CALL HFROTMIX(FMAT, OFMAT, NSPINS, NSBASIS, F, R1, R2, WORK)
        ELSE
            CALL HFLINMIX(FMAT, OFMAT, NSPINS, NSBASIS, F, R1, R2, WORK)
        END IF
        IHFIT = IHFIT + 1
        IF (IHFIT > NHFIT) THEN
            WRITE(stdout, *) "*** WARNING Hartree-Fock did not converge ***"
            BR = .FALSE.
        END IF
        IF (ABS(ECUR - ELAST) <= EDELTA .AND. RMSD <= CDELTA .AND. IHFIT > 5) THEN
            WRITE(stdout, *) "*** Hartree-Fock converged in ", IHFIT, " iterations."
            WRITE(stdout, *) "*** HF ENERGY=", ECUR
            BR = .FALSE.
        END IF
        ELAST = ECUR

        END DO
        IF (BTEST(ILOGGING, 11)) then
            CALL WRITEHFPSIALL(NBASIS, FMAT, HFES, G1, NSPINS, NSBASIS, .FALSE.)
        end if
!.. We write out HFMAT
        HFBASIS = (0.0_dp)
        DO I = 1, NSBASIS
        DO ISPN = 1, NSPINS
            K = (I - 1) * NSPINS + ISPN
            DO J = 1, NSBASIS
                L = (J - 1) * NSPINS + ISPN
!.. eigenvector N is in FMAT(i,N,ISPN), where i is the component of the
!.. vector
!.. HFBASIS(HFBASISFN,PRIMBASISFN) has PRIMBASISFN varying slowest
                HFBASIS(K, L) = FMAT(J, I, ISPN)
            END DO
            HFE(K) = HFES(I, ISPN)
        END DO
        END DO
    END subroutine UHFSCF