subroutine CalcFOCKMatrix()
use SystemData, only: nBasis
use LoggingData, only: tRDMonfly
integer :: i, j, k, l, a, b, ierr
real(dp) :: FOCKDiagSumHF, FOCKDiagSumNew
character(len=*), parameter :: this_routine = 'CalcFOCKMatrix'
!NEED TO FIX THIS!
! This subroutine calculates and writes out the fock matrix for the transformed orbitals.
! ARR is originally the fock matrix in the HF basis.
! ARR(:,1) - ordered by energy, ARR(:,2) - ordered by spin-orbital index.
! When transforming the orbitals into approximate natural orbitals, we want to save memory, so don't bother
! calculating the whole matrix, just the diagonal elements that we actually need.
if (tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) then
allocate(ArrDiagNew(NoOrbs), stat=ierr)
call LogMemAlloc('ArrDiagNew', NoOrbs, 8, this_routine, ArrDiagNewTag, ierr)
ArrDiagNew(:) = 0.0_dp
else
allocate(ArrNew(NoOrbs, NoOrbs), stat=ierr)
call LogMemAlloc('ArrNew', NoOrbs**2, 8, this_routine, ArrNewTag, ierr)
ArrNew(:, :) = 0.0_dp
end if
! First calculate the sum of the diagonal elements, ARR.
! Check if this is already being done.
FOCKDiagSumHF = 0.0_dp
do a = 1, nBasis
FOCKDiagSumHF = FOCKDiagSumHF + Arr(a, 2)
end do
write(stdout, *) 'Sum of the fock matrix diagonal elements in the HF basis set = ', FOCKDiagSumHF
FOCKDiagSumNew = 0.0_dp
do j = 1, NoRotOrbs
l = SymLabelList3_rot(j)
if (tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) then
do a = 1, NoOrbs
b = SymLabelList2_rot(a)
if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
ArrDiagNew(l) = ArrDiagNew(l) + (CoeffT1(a, j) * ARR(b, 2) * CoeffT1(a, j))
else
ArrDiagNew(l) = ArrDiagNew(l) + (CoeffT1(a, j) * ARR(2 * b, 2) * CoeffT1(a, j))
end if
end do
if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew(l))
else
FOCKDiagSumNew = FOCKDiagSumNew + (ArrDiagNew(l) * 2)
end if
else
do i = 1, NoRotOrbs
k = SymLabelList2_rot(i)
ArrNew(k, l) = 0.0_dp
do a = 1, NoOrbs
b = SymLabelList2_rot(a)
if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
ArrNew(k, l) = ArrNew(k, l) + (CoeffT1(a, i) * Arr(b, 2) * CoeffT1(a, j))
else
ArrNew(k, l) = ArrNew(k, l) + (CoeffT1(a, i) * Arr(2 * b, 2) * CoeffT1(a, j))
end if
end do
end do
if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
FOCKDiagSumNew = FOCKDiagSumNew + (ArrNew(l, l))
else
FOCKDiagSumNew = FOCKDiagSumNew + (ArrNew(l, l) * 2)
end if
! Only running through spat orbitals, count each twice to
! compare to above.
end if
end do
! If we are truncation the virtual space, only the unfrozen entries
! will be transformed.
write(stdout, *) 'Sum of the fock matrix diagonal elements in the transformed basis set = ', FOCKDiagSumNew
! Refill ARR(:,1) (ordered in terms of energies), and ARR(:,2) (ordered in terms of orbital number).
! ARR(:,2) needs to be ordered in terms of symmetry and then energy (like SymLabelList), so currently this ordering will not be
! correct when reading in qchem intDUMPS as the orbital number ordering is by energy.
! If we are only writing out 1 ROFCIDUMP or we are not truncating at
! all - can refill ARR etc.
if (NoDumpTruncs <= 1) then
if (tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) then
if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
do j = 1, NoOrbs
ARR(j, 2) = ArrDiagNew(j)
ARR(j, 1) = ArrDiagNew(BRR(j))
end do
else
do j = 1, NoOrbs
ARR(2 * j, 2) = ArrDiagNew(j)
ARR(2 * j - 1, 2) = ArrDiagNew(j)
ARR(2 * j, 1) = ArrDiagNew(BRR(2 * j) / 2)
ARR(2 * j - 1, 1) = ArrDiagNew(BRR(2 * j) / 2)
end do
end if
else
if (tStoreSpinOrbs .or. tTurnStoreSpinOff) then
do j = 1, NoRotOrbs
ARR(j, 2) = ArrNew(j, j)
ARR(j, 1) = ArrNew(BRR(j), BRR(j))
end do
else
do j = 1, NoRotOrbs
ARR(2 * j, 2) = ArrNew(j, j)
ARR(2 * j - 1, 2) = ArrNew(j, j)
ARR(2 * j, 1) = ArrNew(BRR(2 * j) / 2, BRR(2 * j) / 2)
ARR(2 * j - 1, 1) = ArrNew(BRR(2 * j) / 2, BRR(2 * j) / 2)
end do
end if
end if
end if
if ((tUseMP2VarDenMat .or. tFindCINatOrbs .or. tUseHFOrbs .or. tRDMonfly) .and. (NoDumpTruncs <= 1)) then
deallocate(ArrDiagNew)
call LogMemDealloc(this_routine, ArrDiagNewTag)
else if (NoDumpTruncs <= 1) then
deallocate(ArrNew)
call LogMemDealloc(this_routine, ArrNewTag)
end if
write(stdout, *) 'end of calcfockmatrix'
call neci_flush(stdout)
end subroutine CalcFOCKMatrix