subroutine SetUpSymLabels_RDM()
! This routine just sets up the symmetry labels so that the orbitals
! are ordered according to symmetry (all beta then all alpha if spin orbs).
use rdm_data, only: tOpenShell
use RotateOrbsData, only: SymLabelList2_rot, SymLabelCounts2_rot, SymLabelListInv_rot
use RotateOrbsData, only: NoOrbs, SpatOrbs, NoSymLabelCounts
use sort_mod, only: sort
use SystemData, only: G1, BRR, lNoSymmetry, tFixLz, iMaxLz
use UMatCache, only: gtID
use MemoryManager, only: LogMemAlloc, LogMemDealloc
integer, allocatable :: SymOrbs_rot(:)
integer :: SymOrbs_rotTag, ierr, i, j, SpatSym, LzSym
integer :: lo, hi, Symi, SymCurr, Symi2, SymCurr2
character(len=*), parameter :: t_r = 'SetUpSymLabels_RDM'
! This is only allocated temporarily to be used to order the orbitals by.
allocate(SymOrbs_rot(NoOrbs), stat=ierr)
call LogMemAlloc('SymOrbs_rot', NoOrbs, 4, t_r, SymOrbs_rotTag, ierr)
if (ierr /= 0) call stop_all(t_r, "Mem allocation for SymOrbs_rot failed.")
! Now we want to put the spatial orbital index, followed by the symmetry.
SymLabelList2_rot(:) = 0
SymOrbs_rot(:) = 0
! *** STEP 1 *** Fill SymLabelList2_rot.
! Find the orbitals and order them in terms of symmetry.
do i = 1, SpatOrbs
if (tOpenShell) then
! For open shell systems, all alpha are followed by all beta.
SymLabelList2_rot(i) = BRR(2 * i)
SymLabelList2_rot(i + SpatOrbs) = BRR((2 * i) - 1)
if (tFixLz) then
SpatSym = int(G1(BRR(2 * i))%sym%S)
LzSym = int(G1(BRR(2 * i))%Ml)
SymOrbs_rot(i) = (SpatSym * ((2 * iMaxLz) + 1)) + (LzSym + iMaxLz)
SpatSym = int(G1(BRR((2 * i) - 1))%sym%S)
LzSym = int(G1(BRR((2 * i) - 1))%Ml)
SymOrbs_rot(i + SpatOrbs) = (SpatSym * ((2 * iMaxLz) + 1)) + (LzSym + iMaxLz)
else
SymOrbs_rot(i) = int(G1(BRR(2 * i))%sym%S)
SymOrbs_rot(i + SpatOrbs) = int(G1(BRR((2 * i) - 1))%sym%S)
end if
else
SymLabelList2_rot(i) = gtID(BRR(2 * i))
if (tFixLz) then
SpatSym = int(G1(BRR(2 * i))%sym%S)
LzSym = int(G1(BRR(2 * i))%Ml)
SymOrbs_rot(i) = (SpatSym * ((2 * iMaxLz) + 1)) + (LzSym + iMaxLz)
else
SymOrbs_rot(i) = int(G1(BRR(2 * i))%sym%S)
end if
! Orbital BRR(2*i) for i = 1 will be the beta orbital with the
! second lowest energy - want the spatial orbital index to go with this.
! G1 also in spin orbitals - get symmetry of this beta orbital, will
! be the same as the spatial orbital.
end if
end do
call sort(SymOrbs_rot(1:SpatOrbs), SymLabelList2_rot(1:SpatOrbs))
! Sorts SymLabelList2_rot according to the order of SymOrbs_rot
! (i.e. in terms of symmetry).
if (tOpenShell) &
call sort(SymOrbs_rot(SpatOrbs + 1:nBasis), SymLabelList2_rot(SpatOrbs + 1:nBasis))
! Also do this for the beta set if spin orbitals.
! *** STEP 2 *** Fill SymLabelCounts2_rot_rot. This is like
! SymLabelCounts2_rot, but has a different ordering - BEWARE.
! SymLabelCounts(1,:) contains the position in SymLabelList2_rot
! where the symmetry index starts, SymLabelCounts(2,:) contains the
! number of orbitals in that symmetry index. Again if spin orbs, all
! alpha are followed by all beta - i.e. first 8 refer to alpha, second
! 8 to beta.
if (lNoSymmetry) then
! If we are ignoring symmetry, all orbitals essentially have
! symmetry 0.
SymLabelCounts2_rot(1, 1) = 1
SymLabelCounts2_rot(2, 1) = SpatOrbs
if (tOpenShell) then
SymLabelCounts2_rot(1, 9) = SpatOrbs + 1
SymLabelCounts2_rot(2, 9) = SpatOrbs
end if
else
! Otherwise we run through the occupied orbitals, counting the
! number with each symmetry (spatial and Lz) and noting where in
! SymLabelList2_rot each symmetry block starts.
SymCurr = 0
SymLabelCounts2_rot(1, 1) = 1
if (tOpenShell) then
SymCurr2 = 0
SymLabelCounts2_rot(1, 9) = SpatOrbs + 1
end if
do i = 1, SpatOrbs
if (tOpenShell) then
Symi = SymOrbs_rot(i)
Symi2 = SymOrbs_rot(i + SpatOrbs)
else
Symi = SymOrbs_rot(i)
end if
SymLabelCounts2_rot(2, (Symi + 1)) = SymLabelCounts2_rot(2, (Symi + 1)) + 1
if (Symi /= SymCurr) then
do j = SymCurr + 1, Symi
SymLabelCounts2_rot(1, (j + 1)) = i
end do
SymCurr = Symi
end if
if (tOpenShell) then
SymLabelCounts2_rot(2, (Symi2 + 9)) = SymLabelCounts2_rot(2, (Symi2 + 9)) + 1
if (Symi2 /= SymCurr2) then
do j = SymCurr2 + 1, Symi2
SymLabelCounts2_rot(1, (j + 9)) = i + SpatOrbs
end do
SymCurr2 = Symi2
end if
end if
end do
end if
! Go through each symmetry group, making sure the orbitals are
! ordered lowest to highest within each symmetry.
do i = 1, NoSymLabelCounts
if (SymLabelCounts2_rot(2, i) /= 0) then
lo = SymLabelCounts2_rot(1, i)
hi = lo + SymLabelCounts2_rot(2, i) - 1
call sort(SymLabelList2_rot(lo:hi))
end if
end do
! Construct the inverse matrix. While SymLabelList2_rot takes a
! position and tells us what orbital is in it, we also might need to
! take an orbital and find out what position to put its contribution in.
do i = 1, NoOrbs
SymLabelListInv_rot(SymLabelList2_rot(i)) = i
end do
! Deallocate the arrays just used in this routine.
deallocate(SymOrbs_rot)
call LogMemDealloc(t_r, SymOrbs_rotTag)
end subroutine SetUpSymLabels_RDM