subroutine DiagRDM(rdm, SumTrace)
! 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 rdm_data, only: tOpenShell, one_rdm_t
use MemoryManager, only: LogMemAlloc, LogMemDealloc
use SystemData, only: G1, tUseMP2VarDenMat, tFixLz, iMaxLz
type(one_rdm_t), intent(inout) :: rdm
real(dp), intent(out) :: SumTrace
real(dp) :: SumDiagTrace
real(dp), allocatable :: WORK2(:), EvaluesSym(:), NOMSym(:, :)
integer :: ierr, i, j, spin, Sym, LWORK2, WORK2Tag, SymStartInd, NoSymBlock
integer :: EvaluesSymTag, NOMSymTag, k, MaxSym
logical :: tDiffSym, tDiffLzSym
character(len=*), parameter :: t_r = 'DiagRDM'
! Test that we're not breaking symmetry.
! And calculate the trace at the same time.
SumTrace = 0.0_dp
do i = 1, NoOrbs
do j = 1, NoOrbs
tDiffSym = .false.
tDiffLzSym = .false.
if (tOpenShell) then
if ((int(G1(SymLabelList2_rot(i))%sym%S) /= &
int(G1(SymLabelList2_rot(j))%sym%S))) tDiffSym = .true.
if ((int(G1(SymLabelList2_rot(i))%Ml) /= &
int(G1(SymLabelList2_rot(j))%Ml))) tDiffLzSym = .true.
else
if ((int(G1(2 * SymLabelList2_rot(i))%sym%S) /= &
int(G1(2 * SymLabelList2_rot(j))%sym%S))) tDiffSym = .true.
if ((int(G1(2 * SymLabelList2_rot(i))%Ml) /= &
int(G1(2 * SymLabelList2_rot(j))%Ml))) tDiffLzSym = .true.
end if
if (tDiffSym) then
if (abs(rdm%matrix(i, j)) >= 1.0E-15_dp) then
write(stdout, '(6A8,A20)') 'i', 'j', 'Label i', 'Label j', 'Sym i', &
'Sym j', 'Matrix value'
if (tOpenShell) then
write(stdout, '(6I3,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
int(G1(SymLabelList2_rot(i))%sym%S), &
int(G1(SymLabelList2_rot(j))%sym%S), rdm%matrix(i, j)
else
write(stdout, '(6I3,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
int(G1(2 * SymLabelList2_rot(i))%sym%S), &
int(G1(2 * SymLabelList2_rot(j))%sym%S), rdm%matrix(i, j)
end if
if (tUseMP2VarDenMat) then
write(stdout, *) '**WARNING** - There is a non-zero 1-RDM &
&value between orbitals of different symmetry.'
write(stdout, *) 'These elements will be ignored, and the symmetry &
&maintained in the final transformation matrix.'
else
write(stdout, *) 'k,SymLabelList2_rot(k),SymLabelListInv_rot(k)'
do k = 1, NoOrbs
write(stdout, *) k, SymLabelList2_rot(k), SymLabelListInv_rot(k)
end do
call neci_flush(stdout)
call Stop_All(t_r, 'Non-zero rdm%matrix value between &
&different symmetries.')
end if
end if
rdm%matrix(i, j) = 0.0_dp
end if
if (tDiffLzSym) then
if (abs(rdm%matrix(i, j)) >= 1.0E-15_dp) then
write(stdout, '(6A8,A40)') 'i', 'j', 'Label i', 'Label j', 'Lz i', &
'Lz j', 'Matrix value'
if (tOpenShell) then
write(stdout, '(6I8,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
int(G1(SymLabelList2_rot(i))%Ml), &
int(G1(SymLabelList2_rot(j))%Ml), rdm%matrix(i, j)
else
write(stdout, '(6I8,F40.20)') i, j, SymLabelList2_rot(i), SymLabelList2_rot(j), &
int(G1(2 * SymLabelList2_rot(i))%Ml), &
int(G1(2 * SymLabelList2_rot(j))%Ml), rdm%matrix(i, j)
end if
write(stdout, '(A)') ' **WARNING** - There is a non-zero 1-RDM element &
&between orbitals of different Lz symmetry.'
end if
end if
end do
SumTrace = SumTrace + rdm%matrix(i, i)
end do
write(stdout, *) ''
write(stdout, *) 'Calculating eigenvectors and eigenvalues of the 1-RDM'
call neci_flush(stdout)
! 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.
Sym = 0
LWORK2 = -1
if (tOpenShell) then
if (tFixLz) then
MaxSym = (16 * ((2 * iMaxLz) + 1)) - 1
else
MaxSym = 15
end if
else
if (tFixLz) then
MaxSym = (8 * ((2 * iMaxLz) + 1)) - 1
else
MaxSym = 7
end if
end if
do while (Sym <= MaxSym)
NoSymBlock = SymLabelCounts2_rot(2, Sym + 1)
SymStartInd = SymLabelCounts2_rot(1, Sym + 1) - 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(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) = rdm%matrix(SymStartInd + i, SymStartInd + j)
end do
end do
call dsyev('V', 'L', NoSymBlock, NOMSym, NoSymBlock, EvaluesSym, WORK2, LWORK2, ierr)
! NOMSym goes in as the original NOMSym, comes out as the
! eigenvectors (Coefficients).
! EvaluesSym comes out as the eigenvalues in ascending order.
do i = 1, NoSymBlock
rdm%evalues(SymStartInd + i) = EvaluesSym(NoSymBlock - i + 1)
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.
do j = 1, NoSymBlock
do i = 1, NoSymBlock
rdm%matrix(SymStartInd + i, SymStartInd + j) = NOMSym(i, NoSymBlock - j + 1)
end do
end do
! Directly fill the coefficient matrix with the eigenvectors from
! the diagonalization.
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.
rdm%evalues(SymStartInd + 1) = rdm%matrix(SymStartInd + 1, SymStartInd + 1)
rdm%matrix(SymStartInd + 1, SymStartInd + 1) = 1.0_dp
end if
Sym = Sym + 1
end do
write(stdout, *) 'Matrix diagonalised'
call neci_flush(stdout)
SumDiagTrace = 0.0_dp
do i = 1, NoOrbs
SumDiagTrace = SumDiagTrace + rdm%evalues(i)
end do
if ((abs(SumDiagTrace - SumTrace)) > 1.0_dp) then
write(stdout, *) 'Sum of diagonal 1-RDM elements : ', SumTrace
write(stdout, *) 'Sum of eigenvalues : ', SumDiagTrace
write(stdout, *) 'WARNING : &
&The trace of the 1RDM matrix before diagonalisation is '
write(stdout, *) 'not equal to that after.'
end if
! The MO's still correspond to SymLabelList2_rot.
! Although the NO's are slightly jumbled, they are only jumbled within
! their symmetry blocks. They still correspond to the symmetries of
! SymLabelList2_rot, which is the important part.
! But in order to look at the output, it is easier to consider them in
! terms of highest occupied to lowest occupied - i.e. in terms of the
! NO eigenvalues (occupation numbers).
call order_one_rdm(rdm)
end subroutine DiagRDM