subroutine Diagonalizehij()
! This routine takes the original <i|h|j> matrix and diagonalises it. The resulting coefficients from this process
! are then the rotation coefficients to be applied to the four index integrals etc.
! This eliminates the <i|h|j> elements from the single excitations, and leaves only coulomb and exchange terms.
! In order to maintain the same HF energy, only the virtual elements are diagonalised, within symmetry blocks.
integer :: i, j, Sym, ierr, NoSymBlock, WorkSize, WorkCheck, SymStartInd
integer(TagIntType) WorkTag, DiagTMAT2DBlockTag, TMAT2DSymBlockTag
real(dp), allocatable :: TMAT2DSymBlock(:, :), DiagTMAT2DBlock(:), Work(:)
character(len=*), parameter :: this_routine = 'Diagonalizehij'
write(stdout, *) 'The original coefficient matrix'
do i = 1, NoOrbs
do j = 1, NoOrbs
write(stdout, '(F20.10)', advance='no') CoeffT1(j, i)
end do
write(stdout, *) ''
end do
write(stdout, *) 'The original TMAT2D matrix'
do i = 1, NoOrbs
do j = 1, NoOrbs
write(stdout, '(F20.10)', advance='no') TMAT2DTemp(j, i)
end do
write(stdout, *) ''
end do
TMAT2DRot(:, :) = 0.0_dp
DiagTMAT2Dfull(:) = 0.0_dp
! Now need to pick out symmetry blocks, from the virtual orbitals and diagonalize them.
! Take first symmetry, (0) and find the number of virtual orbitals with this symmetry. If this is greater than 1,
! take the block, diagonlize it, and put it into TMAT2DRot.
Sym = 0
WorkSize = -1
do while (Sym <= 7)
NoSymBlock = SymLabelCounts2_rot(2, Sym + 9)
SymStartInd = SymLabelCounts2_rot(1, Sym + 9) - 1
! This is one less than the index that the symmetry starts, so that when we run through i = 1,..., we can
! start at SymStartInd+i.
if (NoSymBlock > 1) then
allocate(TMAT2DSymBlock(NoSymBlock, NoSymBlock), stat=ierr)
call LogMemAlloc('TMAT2DSymBlock', NoSymBlock**2, 8, this_routine, TMAT2DSymBlockTag, ierr)
allocate(DiagTMAT2DBlock(NoSymBlock), stat=ierr)
call LogMemAlloc('DiagTMAT2DBlock', NoSymBlock, 8, this_routine, DiagTMAT2DBlockTag, ierr)
WorkCheck = 3 * NoSymBlock + 1
WorkSize = WorkCheck
allocate(Work(WorkSize), stat=ierr)
call LogMemAlloc('Work', WorkSize, 8, this_routine, WorkTag, ierr)
do j = 1, NoSymBlock
do i = 1, NoSymBlock
TMAT2DSymBlock(i, j) = TMAT2DTemp(SymStartInd + i, SymStartInd + j)
end do
end do
write(stdout, *) '*****'
write(stdout, *) 'Symmetry ', Sym, ' has ', NoSymBlock, ' orbitals .'
write(stdout, *) 'The TMAT2D for this symmetry block is '
do i = 1, NoSymBlock
do j = 1, NoSymBlock
write(stdout, '(F20.10)', advance='no') TMAT2DSymBlock(j, i)
end do
write(stdout, *) ''
end do
call DSYEV('V', 'U', NoSymBlock, TMAT2DSymBlock, NoSymBlock, DiagTMAT2Dblock, Work, WorkSize, ierr)
! TMAT2DSymBlock goes in as the original TMAT2DSymBlock, comes out as the eigenvectors (Coefficients).
! TMAT2DBlock comes out as the eigenvalues in ascending order.
if (ierr /= 0) then
write(stdout, *) 'Problem with symmetry, ', Sym, ' of TMAT2D'
call neci_flush(stdout)
call Stop_All(this_routine, "Diagonalization of TMAT2DSymBlock failed...")
end if
write(stdout, *) 'After diagonalization, the e-vectors (diagonal elements) of this matrix are,'
do i = 1, NoSymBlock
write(stdout, '(F20.10)', advance='no') DiagTMAT2Dblock(i)
end do
write(stdout, *) ''
write(stdout, *) 'These go from orbital,', SymStartInd + 1, ' to ', SymStartInd + NoSymBlock
do i = 1, NoSymBlock
DiagTMAT2Dfull(SymStartInd + i - NoOcc) = DiagTMAT2DBlock(i)
end do
! CAREFUL if eigenvalues are put in ascending order, this may not be correct, with the labelling system.
! may be better to just take coefficients and transform TMAT2DRot in transform2elints.
! a check that comes out as diagonal is a check of this routine anyway.
write(stdout, *) 'The eigenvectors (coefficients) for symmtry block ', Sym
do i = 1, NoSymBlock
do j = 1, NoSymBlock
write(stdout, '(F20.10)', advance='no') TMAT2DSymBlock(j, i)
end do
write(stdout, *) ''
end do
! Directly fill the coefficient matrix with the eigenvectors from the diagonalization.
do j = 1, NoSymBlock
do i = 1, NoSymBlock
CoeffT1(SymStartInd + i, SymStartInd + j) = TMAT2DSymBlock(i, j)
end do
end do
deallocate(Work)
call LogMemDealloc(this_routine, WorkTag)
deallocate(DiagTMAT2DBlock)
call LogMemDealloc(this_routine, DiagTMAT2DBlockTag)
deallocate(TMAT2DSymBlock)
call LogMemDealloc(this_routine, TMAT2DSymBlockTag)
else if (NoSymBlock == 1) then
DiagTMAT2Dfull(SymStartInd + 1 - NoOcc) = TMAT2DTemp(SymStartInd + 1, SymStartInd + 1)
write(stdout, *) '*****'
write(stdout, *) 'Symmetry ', Sym, ' has only one orbital.'
write(stdout, *) 'Copying diagonal element,', SymStartInd + 1, 'to DiagTMAT2Dfull'
end if
Sym = Sym + 1
end do
write(stdout, *) '*****'
write(stdout, *) 'The final coefficient matrix'
do i = 1, NoOrbs
do j = 1, NoOrbs
write(stdout, '(F20.10)', advance='no') CoeffT1(j, i)
end do
write(stdout, *) ''
end do
write(stdout, *) '*****'
write(stdout, *) 'The diagonal elements of TMAT2D'
do i = 1, (NoOrbs - NoOcc)
write(stdout, *) DiagTMAT2Dfull(i)
end do
end subroutine Diagonalizehij