SUBROUTINE UHFGRADDESC(NBASIS, NBASISMAX, G1, BRR, ECORE, UMAT, HFE, HFBASIS, &
NHFIT, NEL, MS, NSPINS, NSBASIS, HFES, HFMIX, CMAT, OCMAT, DEDCIJ, &
DMAT, EDELTA, CDELTA, R1, R2, WORK, TRHF, IHFMETHOD, &
TREADHF, FRAND, HFDET, ILOGGING)
TYPE(BasisFn) G1(*)
INTEGER NSPINS, NSBASIS
INTEGER NBASIS, nBasisMax(5, *)
real(dp) UMAT(*)
real(dp) ECORE
INTEGER BRR(NBASIS)
real(dp) HFBASIS(NBASIS, NBASIS), HFE(NBASIS)
real(dp) HFES(NSBASIS, NSPINS)
HElement_t(dp) CMAT(NSBASIS, NSBASIS, NSPINS)
real(dp) OCMAT(NSBASIS, NSBASIS, NSPINS)
HElement_t(dp) DEDCIJ(NSBASIS, NSBASIS, NSPINS)
HElement_t(dp) DMAT(NSBASIS, NSBASIS, NSPINS)
HElement_t(dp) WORK(NBASIS * 3)
INTEGER INORDER(100, 2), ILOGGING
real(dp) EORDER(100, 2)
!,HFMIX
HElement_t(dp) R1(NSBASIS, NSBASIS), R2(NSBASIS, NSBASIS)
INTEGER NHFIT, NEL, IHFMETHOD, NELEX2
!,NSTART(NEL)
INTEGER I, J, K, L, ISPN, NELS(NSPINS)
real(dp) ELAST, HFMIX, ECUR, FRAND
! INTEGER INDS(NBASIS)
! real(dp) TOT,ELAST,EDELTA,ECUR,TOT2,CDELTA
! INTEGER ID1,ID2,ID3,ID4,INFO,
INTEGER IHFIT
INTEGER MS
! real(dp) F
! real(dp) RMSD
LOGICAL BR
LOGICAL TRHF, TREADHF
real(dp) TOT, MIX, TOT2
INTEGER NDET1(0:NEL + 1), NSPN(NSPINS), NELEX
INTEGER IRHFB, JSPN
real(dp) RMSD, EDELTA, CDELTA, EN, ECUR2
INTEGER HFDET(*)
character(*), parameter :: this_routine = 'UHFGRADDESC'
irhfb = 0
IF (NSBASIS > 99) call stop_all(this_routine, 'ERROR - hardcoded NSBASIS limit of 100')
ELAST = 1.D20
WRITE(stdout, *) "Performing Hartree-Fock Gradient Descent..."
IF (IHFMETHOD == 1) THEN
WRITE(stdout, *) "Method 1:Singles replacement "
ELSEIF (IHFMETHOD == 2) THEN
WRITE(stdout, *) "Method 2:Explicit differential"
END IF
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
!.. Cij - i corresponds to rows and new basis functions, phi_i
!.. j corresponds to columns and old basis functions, u_j
!.. phi_i=sum_j=1,M cij u_j
! VMAT=0.0_dp
IF (TREADHF) THEN
CALL READHFFMAT(NBASIS, CMAT, HFES, G1, NSPINS, NSBASIS, .TRUE.)
ELSE
CALL GENHFGUESS(CMAT, NSPINS, NSBASIS, BRR, G1, .TRUE., MS, FRAND, NELS, HFDET)
END IF
! DO ISPN=1,NSPINS
! DO I=1,NSBASIS
! WRITE(stdout,*) (FMAT(I,J,ISPN),J=1,NSBASIS)
! ENDDO
! ENDDO
!.. Initialize our HF det in NDET1
DO I = 1, NSPINS
NSPN(I) = 0
END DO
I = 1
NDET1(0) = 0
NDET1(NEL + 1) = NBASIS + 1
DO WHILE (I <= NEL)
DO ISPN = 1, NSPINS
IF (NSPN(ISPN) < NELS(ISPN)) THEN
NDET1(I) = NSPN(ISPN) * NSPINS + ISPN
NSPN(ISPN) = NSPN(ISPN) + 1
I = I + 1
END IF
END DO
END DO
call write_det(stdout, NDET1(1), .true.)
BR = .TRUE.
WRITE(stdout, *) "Iteration Energy MSD Fock Energy"
IHFIT = 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, CMAT(1, 1, IRHFB), 1, CMAT(1, 1, 3 - IRHFB), 1)
END IF
CALL DCOPY(NSBASIS * NSBASIS * NSPINS, CMAT, 1, OCMAT, 1)
IHFIT = IHFIT + 1
!.. First Calculate dE/dcij
!.. dE/dcij is automatically 0 if i>N as the HF det only depends on
!.. phi_1 to phi_N. All values of j must be iterated as each phi_i is
!.. dependent on all u_j
ECUR = GETHELEMENT2T(NDET1(1), NDET1(1), NEL, NBASISMAX, NBASIS, ECORE, 0, CMAT, NSBASIS, NSPINS)
!.. Calculate the Gradient
IF (IHFMETHOD == 1) THEN
CALL CALCDEDCIJ(CMAT, DEDCIJ, NDET1, NSPINS, NSBASIS, ECORE, NBASIS, NBASISMAX, NELS, NEL, ECUR)
ELSEIF (IHFMETHOD == 2) THEN
CALL CALCDEDCIJ2(CMAT, DEDCIJ, NDET1, NSPINS, NSBASIS, UMAT, NELS, NEL)
END IF
!.. DEDCIJ now comtains all elements of dE/dcij
!.. To move down the slope, we subtract a small amount of this from cij,
!.. and re-orthogonalise
!.. HFMIX is -ve
! WRITE(stdout,*) ((CMAT(I,J,ISPN),I=1,NSBASIS),J=1,NSBASIS)
! WRITE(stdout,*)
! WRITE(stdout,*) ((DEDCIJ(I,J,ISPN),I=1,NSBASIS),J=1,NSBASIS)
!.. modify the velocity.
! R1 = VMAT(:,:,ISPN) + MIX*DEDCIJ(:,:,ISPN)
! CALL DCOPY(NSBASIS*NSBASIS,R1,1,VMAT(1,1,ISPN),1)
! R1 = CMAT(:,:,ISPN) + 0.1_dp*VMAT(:,:,ISPN)
!.. Remove the projection of the "force" already in the direction of
!.. the coefficients
DO ISPN = 1, NSPINS
R1 = 0.0_dp
DO I = 1, NELS(ISPN)
TOT = 0.0_dp
DO J = 1, NSBASIS
TOT = TOT + DEDCIJ(I, J, ISPN) * CMAT(I, J, ISPN)
END DO
TOT2 = 0.0_dp
DO J = 1, NSBASIS
R1(I, J) = DEDCIJ(I, J, ISPN) - TOT * CMAT(I, J, ISPN)
TOT2 = TOT2 + R1(I, J)**2
END DO
TOT2 = SQRT(TOT2)
DO J = 1, NSBASIS
! R1(I,J)=R1(I,J)/TOT2
END DO
END DO
MIX = -HFMIX
!/ABS(ECUR)
! IF(ABS(ECUR).LT.1.0e-4_dp) MIX=HFMIX
IF (IHFIT > NHFIT) BR = .FALSE.
R2 = CMAT(:, :, ISPN) + R1
CALL DCOPY(NSBASIS * NSBASIS, R2, 1, CMAT(1, 1, ISPN), 1)
! WRITE(stdout,*)
! WRITE(stdout,*) ((CMAT(I,J,ISPN),I=1,NSBASIS),J=1,NSBASIS)
! WRITE(stdout,*)
! WRITE(stdout,*)
CALL LOWDIN_ORTH(CMAT(1, 1, ISPN), NSBASIS, R1, R2, WORK)
END DO
RMSD = 0.0_dp
DO ISPN = 1, NSPINS
DO I = 1, NSBASIS
DO J = 1, NSBASIS
IF (.NOT. TRHF .OR. (TRHF .AND. ISPN == IRHFB)) then
RMSD = RMSD + (CMAT(I, J, ISPN) - OCMAT(I, J, ISPN))**2
end if
END DO
END DO
END DO
RMSD = SQRT(RMSD / (NSBASIS * NSBASIS * NSPINS))
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
!.. Construct the Density Matrix
CALL GENDMAT(NSPINS, NSBASIS, NELS, CMAT, DMAT, .TRUE.)
!.. Use the Density Matrix to generate the Fock matrix (in DEDCIJ)
CALL GENFMAT(DEDCIJ, DMAT, NSBASIS, NSPINS)
CALL DIAGFMAT(NSPINS, NSBASIS, NELS, DEDCIJ, DMAT, HFES, WORK, ECORE, ECUR2)
WRITE(stdout, "(I6)", advance='no') IHFIT
WRITE(stdout, *) ECUR, RMSD, ECUR2
!.. DEDCIJ now contains HF orbitals, and ECUR the Fock Energy
!.. eigenvector N is in FMAT(i,N,ISPN), where i is the component of the
!.. vector
!.. calculate the orbital energies every time
DO ISPN = 1, NSPINS
DO I = 1, NSBASIS
NELEX = (I - 1) * NSPINS + 1 + ISPN - 1
EN = 0.0_dp
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL GETTRTMATEL(NELEX, NELEX, CMAT, NSBASIS, NSPINS, EN)
#endif
HFES(I, ISPN) = EN
DO JSPN = 1, NSPINS
DO J = 1, NELS(JSPN)
NELEX2 = (J - 1) * NSPINS + 1 + JSPN - 1
IF (NELEX /= NELEX2) THEN
!.. we're not allowed to count the current electron again
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL GETTRUMATEL(NELEX, NELEX2, NELEX, NELEX2, CMAT, NSBASIS, NSPINS, EN)
#endif
! HFES(I,ISPN)=HFES(I,ISPN)+EN
EN = 0.0_dp
IF (ISPN == JSPN) then
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL GETTRUMATEL(NELEX, NELEX2, NELEX2, NELEX, CMAT, NSBASIS, NSPINS, EN)
#endif
endif
! HFES(I,ISPN)=HFES(I,ISPN)-EN
END IF
END DO
END DO
END DO
DO I = 1, NSBASIS
INORDER(I, ISPN) = I
END DO
CALL DCOPY(NSBASIS, HFES(1, ISPN), 1, EORDER(1, ISPN), 1)
call sort(eorder(1:nsBasis, iSpn), inOrder(1:nsBasis, iSpn))
END DO
IF (BR) THEN
!.. Now re-order the orbitals
!.. although we don't actually use this info
DO ISPN = 1, NSPINS
R1 = 0.0_dp
DO I = 1, NSBASIS
!.. If R1(2,1) is occupied, then postmultiplying C by R1 moves
!.. column 2 in C to column 1.
!.. INORDER(1,ISPN) is the old index of the lowest energy orb
R1(I, INORDER(I, ISPN)) = 1
END DO
!.. Work out R2=1.0_dp CMAT*R1+ 0.0_dp*R2
CALL DGEMM('N', 'N', NSBASIS, NSBASIS, NSBASIS, 1.0_dp, R1, NSBASIS, CMAT(1, 1, ISPN), NSBASIS, 0.0_dp, R2, NSBASIS)
! CALL DCOPY(NSBASIS*NSBASIS,R2,1,CMAT(1,1,ISPN),1)
CALL DGEMM('N', 'N', NSBASIS, NSBASIS, NSBASIS, 1.0_dp, R1, NSBASIS, OCMAT(1, 1, ISPN), NSBASIS, 0.0_dp, R2, NSBASIS)
! CALL DCOPY(NSBASIS*NSBASIS,R2,1,OCMAT(1,1,ISPN),1)
END DO
DO I = 1, NSBASIS
DO ISPN = 1, NSPINS
! WRITE(stdout,*) I,ISPN*2-3,HFES(I,ISPN)
! WRITE(stdout,*) I,ISPN*2-3,HFES(INORDER(I,ISPN),ISPN)
!,
! & INORDER(I,ISPN)
END DO
END DO
DO ISPN = 1, NSPINS
DO I = 1, NSBASIS
! WRITE(stdout,*) I,ISPN*2-3,EORDER(I,ISPN),INORDER(I,ISPN)
!HFES(INORDER(I,ISPN),ISPN),
! & INORDER(I,ISPN)
END DO
END DO
END IF
IF (MOD(IHFIT, 10) == 0 .AND. BTEST(ILOGGING, 11)) then
CALL WRITEHFPSIALL(NBASIS, CMAT, HFES, G1, NSPINS, NSBASIS, .TRUE.)
end if
END DO
IF (BTEST(ILOGGING, 11)) then
CALL WRITEHFPSIALL(NBASIS, CMAT, HFES, G1, NSPINS, NSBASIS, .TRUE.)
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) = CMAT(I, J, ISPN)
END DO
HFE(K) = HFES(I, ISPN)
END DO
END DO
RETURN
END subroutine UHFGRADDESC