subroutine DiagNatOrbMat()
! The diagonalisation routine reorders the orbitals in such a way that
! the corresponding orbital labels are lost. In order to keep the spin
! and spatial symmetries, each symmetry must be fed into the
! diagonalisation routine separately. The best way to do this is to
! order the orbitals so that all the alpha orbitals follow all the beta
! orbitals, with the occupied orbitals first, in terms of symmetry, and
! the virtual second, also ordered by symmetry. This gives us
! flexibility w.r.t rotating only the occupied or only virtual and
! looking at high spin states.
use MemoryManager, only: TagIntType
real(dp) :: SumTrace, SumDiagTrace
real(dp), allocatable :: WORK2(:), EvaluesSym(:), NOMSym(:, :)
integer :: ierr, i, j, x, z, Sym, LWORK2
integer :: SymStartInd, NoSymBlock, PrevSym, StartOccVirt, EndOccVirt, Prev, NoOcc
integer(TagIntType) :: EvaluesSymTag, NOMSymTag, WORK2Tag
character(len=*), parameter :: t_r = 'DiagNatOrbMat'
DiagNatOrbMat_Time%timer_name = 'DiagNatOrb'
call set_timer(DiagNatOrbMat_Time, 30)
do x = 1, NoSpinCyc
if (tSeparateOccVirt) then
if (x == 1) then
if (tStoreSpinOrbs) then
NoOcc = nOccBeta
else
NoOcc = NEl / 2
end if
Prev = 0
else if (x == 2) then
NoOcc = nOccAlpha
Prev = SpatOrbs
end if
else
NoOcc = 0
end if
if (tRotateVirtOnly) then
do i = 1, NoOcc
do j = 1, SpatOrbs
NatOrbMat(i + Prev, j + Prev) = 0.0_dp
NatOrbMat(j + Prev, i + Prev) = 0.0_dp
if (i == j) NatOrbMat(i + Prev, j + Prev) = 1.0_dp
end do
Evalues(i + Prev) = 1.0_dp
end do
else if (tRotateOccOnly) then
do i = NoOcc + 1, SpatOrbs
do j = 1, SpatOrbs
NatOrbMat(i + Prev, j + Prev) = 0.0_dp
NatOrbMat(j + Prev, i + Prev) = 0.0_dp
if (i == j) NatOrbMat(i + Prev, j + Prev) = 1.0_dp
end do
Evalues(i + Prev) = 1.0_dp
end do
else if (tSeparateOccVirt) then
do i = 1, NoOcc
do j = NoOcc + 1, SpatOrbs
NatOrbMat(i + Prev, j + Prev) = 0.0_dp
NatOrbMat(j + Prev, i + Prev) = 0.0_dp
end do
end do
end if
end do
! Test that we're not breaking symmetry.
do i = 1, NoOrbs
do j = 1, NoOrbs
if (tStoreSpinOrbs) then
if ((int(G1(SymLabelList2_rot(i))%sym%S, 4) /= int(G1(SymLabelList2_rot(j))%sym%S, 4))) then
if (abs(NatOrbMat(i, j)) >= 1.0e-15_dp) then
write(stdout, '(6A8,A20)') 'i', 'j', 'Label i', 'Label j', 'Sym i', 'Sym j', 'Matrix value'
write(stdout, '(6I3,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
int(G1(SymLabelList2_rot(i))%sym%S, 4), &
int(G1(SymLabelList2_rot(j))%sym%S, 4), NatOrbMat(i, j)
if (tUseMP2VarDenMat) then
write(stdout, *) '**WARNING** - There is a non-zero NatOrbMat value between " &
& //"orbitals of different symmetry.'
write(stdout, *) 'These elements will be ignored, and the symmetry maintained " &
& //"in the final transformation matrix.'
else
call stop_all(t_r, 'Non-zero NatOrbMat value between different symmetries.')
end if
end if
NatOrbMat(i, j) = 0.0_dp
end if
else
if ((int(G1(SymLabelList2_rot(i) * 2)%sym%S, 4) /= int(G1(SymLabelList2_rot(j) * 2)%sym%S, 4))) then
if (abs(NatOrbMat(i, j)) >= 1.0e-15_dp) then
write(stdout, '(6A8,A20)') 'i', 'j', 'Label i', 'Label j', 'Sym i', 'Sym j', 'Matrix value'
write(stdout, '(6I3,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
int(G1(SymLabelList2_rot(i) * 2)%sym%S, 4), int(G1(SymLabelList2_rot(j) * 2)%sym%S, 4), NatOrbMat(i, j)
if (tUseMP2VarDenMat) then
write(stdout, *) '**WARNING** - There is a non-zero NatOrbMat value between orbitals of " &
& //"different symmetry.'
write(stdout, *) 'These elements will be ignored, and the symmetry maintained in the " &
& //"final transformation matrix.'
else
call stop_all(t_r, 'Non-zero NatOrbMat value between different symmetries.')
end if
end if
NatOrbMat(i, j) = 0.0_dp
end if
end if
end do
end do
SumTrace = 0.0_dp
do i = 1, NoOrbs
SumTrace = SumTrace + NatOrbMat(i, i)
end do
write(stdout, *) 'Calculating eigenvectors and eigenvalues of NatOrbMat'
call neci_flush(stdout)
! If we are using spin orbitals, need to feed in the alpha and beta
! spins separately. Otherwise these jumble up and the final ordering
! is uncorrect. There should be no non-zero values between these, but
! can put a check in for this.
do x = 1, NoSpinCyc
! If we want to maintain the symmetry, we cannot have all the
! orbitals jumbled up when the diagonaliser reorders the
! eigenvectors. Must instead feed each symmetry block in
! separately. This means that although the transformed orbitals
! are jumbled within the symmetry blocks, the symmetry labels are
! all that are relevant and these are unaffected.
StartOccVirt = 1
EndOccVirt = 2
if (tRotateVirtOnly) StartOccVirt = 2
if (tRotateOccOnly) EndOccVirt = 1
do z = StartOccVirt, EndOccVirt
if (x == 1) then
if (z == 1) PrevSym = 1
if (z == 2) PrevSym = 9
else if (x == 2) then
if (z == 1) PrevSym = 17
if (z == 2) PrevSym = 25
end if
Sym = 0
LWORK2 = -1
do while (Sym <= 7)
NoSymBlock = SymLabelCounts2_rot(2, Sym + PrevSym)
! 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.
SymStartInd = SymLabelCounts2_rot(1, Sym + PrevSym) - 1
if (NoSymBlock > 1) then
allocate(NOMSym(NoSymBlock, NoSymBlock), stat=ierr)
call LogMemAlloc('NOMSym', NoSymBlock**2, 8, t_r, NOMSymTag, ierr)
if (ierr /= 0) call stop_all(t_r, "Problem allocating NOMSym.")
allocate(EvaluesSym(NoSymBlock), stat=ierr)
call LogMemAlloc('EvaluesSym', NoSymBlock, 8, t_r, EvaluesSymTag, ierr)
if (ierr /= 0) call stop_all(t_r, "Problem allocating EvaluesSym.")
LWORK2 = 3 * NoSymBlock + 1
allocate(WORK2(LWORK2), stat=ierr)
call LogMemAlloc('WORK2', LWORK2, 8, t_r, WORK2Tag, ierr)
if (ierr /= 0) call stop_all(t_r, "Problem allocating WORK2.")
do j = 1, NoSymBlock
do i = 1, NoSymBlock
NOMSym(i, j) = NatOrbMat(SymStartInd + i, SymStartInd + j)
end do
end do
write(stdout, *) '*****'
write(stdout, *) 'Symmetry ', Sym, 'with x ', x, ' and z ', z, ' has ', NoSymBlock, ' orbitals.'
write(stdout, *) 'The NatOrbMat for this symmetry block is '
do i = 1, NoSymBlock
do j = 1, NoSymBlock
write(stdout, '(F20.10)', advance='no') NOMSym(j, i)
end do
write(stdout, *) ''
end do
! NOMSym goes in as the original NOMSym, comes out as
! the eigenvectors (Coefficients). EvaluesSym comes out
! as the eigenvalues in ascending order.
call dsyev('V', 'L', NoSymBlock, NOMSym, NoSymBlock, EvaluesSym, WORK2, LWORK2, ierr)
write(stdout, *) 'After diagonalization, the e-vectors (diagonal elements) of this matrix are,'
do i = 1, NoSymBlock
write(stdout, '(F20.10)', advance='no') EvaluesSym(i)
end do
write(stdout, *) ''
write(stdout, *) 'These go from orbital,', SymStartInd + 1, ' to ', SymStartInd + NoSymBlock
do i = 1, NoSymBlock
Evalues(SymStartInd + i) = EvaluesSym(i)
end do
! CAREFUL if eigenvalues are put in ascending order,
! this may not be correct, with the labelling system.
! Maybe 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 symmetry block ', Sym
do i = 1, NoSymBlock
do j = 1, NoSymBlock
write(stdout, '(F20.10)', advance='no') NOMSym(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
NatOrbMat(SymStartInd + i, SymStartInd + j) = NOMSym(i, j)
end do
end do
deallocate(WORK2)
call LogMemDealloc(t_r, WORK2Tag)
deallocate(NOMSym)
call LogMemDealloc(t_r, NOMSymTag)
deallocate(EvaluesSym)
call LogMemDealloc(t_r, EvaluesSymTag)
else if (NoSymBlock == 1) then
! The eigenvalue is the lone value, while the
! eigenvector is 1.
Evalues(SymStartInd + 1) = NatOrbMat(SymStartInd + 1, SymStartInd + 1)
NatOrbMat(SymStartInd + 1, SymStartInd + 1) = 1.0_dp
write(stdout, *) '*****'
write(stdout, *) 'Symmetry ', Sym, ' has only one orbital.'
write(stdout, *) 'Copying diagonal element,', SymStartInd + 1, 'to NatOrbMat'
end if
Sym = Sym + 1
end do
end do
end do
write(stdout, *) 'Matrix diagonalised'
call neci_flush(stdout)
SumDiagTrace = 0.0_dp
do i = 1, NoOrbs
SumDiagTrace = SumDiagTrace + Evalues(i)
end do
if ((abs(SumDiagTrace - SumTrace)) > 10.0_dp) then
write(stdout, *) 'Sum of diagonal NatOrbMat elements : ', SumTrace
write(stdout, *) 'Sum of eigenvalues : ', SumDiagTrace
write(stdout, *) 'WARNING, The trace of the 1RDM matrix before diagonalisation is not equal to that after.'
end if
call halt_timer(DiagNatOrbMat_Time)
end subroutine DiagNatOrbMat