subroutine HFDoCalc()
Use global_utilities
use SystemData, only: BasisFN
use IntegralsData, only: tHFBasis, tHFCalc, iHFMethod, tReadHF, nHFIt, HFMix, HFCDelta, HFEDelta
use IntegralsData, only: HFRand, tRHF, ntFrozen, tReadTUMat
use SystemData, only: tCPMD, tHFOrder, nBasisMax, G1, Arr, Brr, ECore, nEl, nBasis, iSpinSkip, LMS
use SystemData, only: tHub, lmsbasis
Use LoggingData, only: iLogging
use IntegralsData, only: UMat, tagUMat, umat_win
Use UMatCache, only: GetUMatSize
Use OneEInts, only: TMat2D, TMat2D2, SetupTMat2, DestroyTMat
use sort_mod
use shared_memory_mpi
use HElem, only: helement_t_size, helement_t_sizeb
use MemoryManager, only: TagIntType
character(25), parameter :: this_routine = 'HFDoCalc'
HElement_t(dp), ALLOCATABLE :: HFBASIS(:), HFE(:)
HElement_t(dp), pointer :: UMat2(:)
INTEGER(MPIArg):: umat2_win
integer i
integer nOrbUsed
integer TMatInt
integer(int64) :: UMatInt
integer(TagIntType), save :: tagUMat2 = 0, tagHFE = 0, tagHFBasis = 0
!C.. If we are using an HF basis instead of our primitive basis, we need
!C.. to load in the coeffs of the HF eigenfunctions in terms of the
!C.. primitive basis.
!C.. We load the coeffs from a file HFBASIS
IF (THFBASIS .OR. THFCALC .OR. (THFORDER .AND. .NOT. TCPMD)) THEN
! NOTE: while HFBasis and HFE are declared to be arrays,
! many of the following routines assume them to be real arrays.
! These need to be changed for use with complex code.
allocate(HFBasis(nBasis * nBasis))
call LogMemAlloc('HFBASIS', nBasis * nBasis, HElement_t_size * 8, this_routine, tagHFBasis)
HFBASIS = (0.0_dp)
!C.. Allocate an array to store the HF Energies
allocate(HFE(nBasis))
call LogMemAlloc('HFE', nBasis, HElement_t_size * 8, this_routine, tagHFE)
HFE = (0.0_dp)
IF (THFORDER .AND. .NOT. THFBASIS) THEN
!C.. If we're not using HF, but just calculating the HF order
!C.. We generate the HF energies (this has no mixing or randomisation, so should jsut
!C.. re-order the orbitals and give us some energy)
!C.. HF basis is NOT using the LMS value set in the input
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL CALCHFBASIS(NBASIS, NBASISMAX, G1, BRR, ECORE, UMAT, HFE, HFBASIS, 1, NEL, LMSBASIS, 1.0_dp, HFEDELTA, &
HFCDELTA, .TRUE., 0, TREADHF, 0.0_dp, FDET, ILOGGING)
CALL ORDERBASISHF(ARR, BRR, HFE, HFBASIS, NBASIS, FDET, NEL)
#endif
else if (THFCALC) THEN
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL CALCHFBASIS(NBASIS, NBASISMAX, G1, BRR, ECORE, UMAT, HFE, HFBASIS, NHFIT, NEL, LMS, HFMIX, HFEDELTA, &
HFCDELTA, TRHF, IHFMETHOD, TREADHF, HFRAND, FDET, ILOGGING)
CALL SETUPHFBASIS(NBASISMAX, G1, NBASIS, HFE, ARR, BRR)
#endif
else if (THFBASIS) THEN
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL READHFBASIS(HFBASIS, HFE, G1, NBASIS)
CALL SETUPHFBASIS(NBASISMAX, G1, NBASIS, HFE, ARR, BRR)
#endif
end if
write(stdout, *) "FINAL HF BASIS"
CALL writebasis(stdout, G1, nBasis, ARR, BRR)
write(stdout, "(A)", advance='no') " Fermi det (D0):"
call write_det(stdout, FDET, .true.)
CALL neci_flush(stdout)
!C.. If in Hubbard, we generate site-spin occupations
IF (THUB) THEN
! Don't think this works
! CALL GENSITESPINOCC(NBASIS,NBASIS/ISPINSKIP,ISPINSKIP,NBASISMAX,G1,NEL,LMS,BRR,HFBASIS)
end if
!C.. We now generate a new U matrix corresponding to the HF basis fns
!C.. This requires a new matrix UMAT2 to be set up, as our HF basis is
!C.. unrestricted
!THIS ROUTINE NO LONGER WORKS WITH NEW TMAT/UMAT MODULARISATION
IF (THFBASIS) THEN
write(stdout, *) "Allocating TMAT2"
CALL SetupTMAT2(nBasis, 2, TMATINT)
NORBUSED = NBASIS - NTFROZEN
IF (TREADTUMAT) THEN
CALL READHFTMAT(NBASIS)
ELSE
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
CALL CALCHFTMAT(NBASIS, HFBASIS, NORBUSED)
#endif
end if
CALL DestroyTMAT(.false.)
TMAT2D => TMAT2D2
NULLIFY (TMAT2D2)
!C.. Allocate the new matrix
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
!C.. We need to pass the TMAT to CALCHFUMAT as TMAT is no longer diagona
!C.. This also modified G1, ARR, BRR
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
IF (TREADTUMAT) THEN
CALL READHFUMAT(UMAT2, NBASIS)
ELSE
CALL CALCHFUMAT(UMAT, UMAT2, NBASIS, HFBASIS, ISPINSKIP, NORBUSED)
end if
#endif
!C.. Now we can remove the old UMATRIX, and set the pointer UMAT to point
!C.. to UMAT2
call shared_sync_mpi(umat2_win)
LogDealloc(tagUMat)
call shared_deallocate_mpi(umat_win, umat)
!Deallocate(UMat)
umat_win = umat2_win
UMat => UMat2
nullify (UMat2)
tagUMat = tagUMat2
tagUMat2 = 0
!C.. spinskip values
ISPINSKIP = 1
NBASISMAX(2, 3) = 1
!C.. Indicate that we're in UHF and that <D0|H|D1>=0 for D1 being a
!C.. single excitation
IF (LMS == 0) THEN
NBASISMAX(4, 5) = 1
DO I = 1, NEL
NUHFDET(I) = BRR(I)
end do
call sort(nUHFDet(1:nel))
ELSE
NBASISMAX(4, 5) = 2
end if
end if
!C.. Now deallocate the HF arrays
deallocate(HFE, HFBasis)
call LogMemDealloc(this_routine, tagHFE)
call LogMemDealloc(this_routine, tagHFBasis)
CALL neci_flush(stdout)
end if
End Subroutine HFDoCalc