SUBROUTINE SetupTMAT(nBASIS, iSS, iSize)
! In:
! nBasis: number of basis functions (orbitals).
! iSS: ratio of nBasis to the number of spatial orbitals.
! iSize: number of elements in TMat/TMatSym.
! Initial allocation of TMat2D or TMatSym (if using symmetry-compressed
! storage of the <i|h|j> integrals).
IMPLICIT NONE
integer Nirrep, nBasis, iSS, nBi, i, basirrep, t, ierr, iState, nStateIrrep
integer iSize, iunit
character(*), parameter :: thisroutine = 'SetupTMAT'
! If this is a CPMD k-point calculation, then we're operating
! under Abelian symmetry: can use George's memory efficient
! TMAT.
if (tCPMD) tCPMDSymTMat = tKP
if (tVASP) tCPMDSymTMat = .true.
IF (tCPMDSymTMat) THEN
! Set up info for indexing scheme (see TMatInd for full description).
Nirrep = NSYMLABELS
nBi = nBasis / iSS
iSize = 0
if (associated(SymLabelIntsCum)) then
deallocate(SymLabelIntsCum)
call LogMemDealloc(thisroutine, tagSymLabelIntsCum)
end if
if (allocated(StateSymMap)) then
deallocate(StateSymMap)
call LogMemDealloc(thisroutine, tagStateSymMap)
end if
if (associated(SymLabelCountsCum)) then
deallocate(SymLabelCountsCum)
call LogMemDealloc(thisroutine, tagSymLabelCountsCum)
end if
allocate(SymLabelIntsCum(nIrrep))
call LogMemAlloc('SymLabelIntsCum', nIrrep, 4, thisroutine, tagSymLabelIntsCum)
allocate(SymLabelCountsCum(nIrrep))
call LogMemAlloc('SymLabelCountsCum', nIrrep, 4, thisroutine, tagSymLabelCountsCum)
allocate(StateSymMap(nBi))
call LogMemAlloc('StateSymMap', nBi, 4, thisroutine, tagStateSymMap)
SYMLABELCOUNTSCUM(1:Nirrep) = 0
SYMLABELINTSCUM(1:Nirrep) = 0
! Ensure no compiler warnings
basirrep = SYMLABELCOUNTS(2, 1)
do i = 1, Nirrep
basirrep = SYMLABELCOUNTS(2, i)
! Block diagonal.
iSize = iSize + (basirrep * (basirrep + 1)) / 2
! # of integrals in this symmetry class and all preceeding
! symmetry classes.
SYMLABELINTSCUM(i) = iSize
! JSS: we no longer need SymLabelCountsCum, but it's a useful test.
! SymLabelCountsCum is the cumulative total # of basis functions
! in the preceeding symmetry classes.
IF (i == 1) THEN
SYMLABELCOUNTSCUM(i) = 0
ELSE
DO t = 1, (i - 1)
SYMLABELCOUNTSCUM(i) = SYMLABELCOUNTSCUM(i) + SYMLABELCOUNTS(2, t)
end do
end if
! JSS: Label states of symmetry i by the order in which they come.
nStateIrrep = 0
do iState = 1, nBi
if (SymClasses(iState) == i) then
nStateIrrep = nStateIrrep + 1
StateSymMap(iState) = nStateIrrep
end if
end do
end do
! The number of basis functions before the last irrep
! plus the number of basis functions in the last irrep
! (which is in basirrep on exiting the above do loop)
! must equal the total # of basis functions.
IF ((SYMLABELCOUNTSCUM(Nirrep) + basirrep) /= nBI) THEN
iunit = get_free_unit()
open(iunit, file="TMATSYMLABEL", status="unknown")
DO i = 1, Nirrep
write(iunit, *) SYMLABELCOUNTSCUM(i)
end do
write(iunit, *) "***************"
write(iunit, *) NBI
CALL neci_flush(iunit)
close(iunit)
call stop_all(thisroutine, 'Not all basis functions found while setting up TMAT')
end if
!iSize=iSize+2
!This is to allow the index of '-1' in the array to give a zero value
!Refer to TMatSym(-1) for when <i|h|j> is zero by symmetry.
allocate(TMATSYM(-1:iSize), STAT=ierr)
Call LogMemAlloc('TMATSym', iSize + 2, HElement_t_size * 8, thisroutine, tagTMATSYM, ierr)
TMATSYM = (0.0_dp)
ELSE
if (tOneElecDiag) then
! In the UEG, the orbitals are eigenfunctions of KE operator, so TMAT is diagonal.
iSize = nBasis
allocate(TMAT2D(nBasis, 1), STAT=ierr)
call LogMemAlloc('TMAT2D', nBasis, HElement_t_size * 8, thisroutine, tagTMat2D)
TMAT2D = (0.0_dp)
else
! Using a square array to hold <i|h|j> (incl. elements which are
! zero by symmetry).
iSize = nBasis * nBasis
allocate(TMAT2D(nBasis, nBasis), STAT=ierr)
call LogMemAlloc('TMAT2D', nBasis * nBasis, HElement_t_size * 8, thisroutine, tagTMat2D)
TMAT2D = (0.0_dp)
end if
end if
END SUBROUTINE SetupTMAT