subroutine EquateDiagFock()
integer :: irr, NumInSym, Orbi, Orbj, w, i, j, k, ConjInd, OrbjConj
real(dp) :: Angle, AngleConj, Check, Norm
CoeffT1(:, :) = 0.0_dp
! Do virtual and occupied orbitals seperately.
do w = MinOccVirt, MaxOccVirt
! Loop over irreps.
do irr = 1, 8
NumInSym = SymLabelCounts2_rot(2, (w - 1) * 8 + irr)
! Loop over the j-orthogonal vectors to create in this
! symmetry block.
do j = 1, NumInSym
Orbj = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + j)
! See if this vector has already been done.
Check = 0.0_dp
do i = 1, NoOrbs
Check = Check + CoeffT1(i, Orbj)
end do
if (.not. near_zero(Check)) then
! This vector is a conjugate pair of another vector and
! has already been worked out...
cycle
end if
! Find out if we this vector will be complex. It will be
! real if j = N or j = N/2
if (j == NumInSym) then
!i The vector will be the normalized 1, 1, 1 vector.
do i = 1, NumInSym
Orbi = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + i)
CoeffT1(Orbi, Orbj) = 1 / SQRT(real(NumInSym, dp))
end do
else if ((mod(NumInSym, 2) == 0) .and. (j == (NumInSym / 2))) then
do i = 1, NumInSym
Orbi = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + i)
if (mod(i, 2) == 1) then
CoeffT1(Orbi, Orbj) = -1 / SQRT(real(NumInSym, dp))
else
CoeffT1(Orbi, Orbj) = 1 / SQRT(real(NumInSym, dp))
end if
end do
else
! Vector is complex - find its conjugate vector - do
! these at the same time.
ConjInd = NumInSym - j
OrbjConj = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + ConjInd)
do i = 1, NumInSym
Orbi = SymLabelList2_rot(SymLabelCounts2_rot(1, (w - 1) * 8 + irr) - 1 + i)
Angle = real(i * j * 2, dp) * PI / real(NumInSym, dp)
AngleConj = real(i * ConjInd * 2, dp) * PI / real(NumInSym, dp)
CoeffT1(Orbi, Orbj) = (1 / SQRT(real(2 * NumInSym, dp))) * (COS(Angle) + COS(AngleConj))
CoeffT1(Orbi, OrbjConj) = (1 / SQRT(real(2 * NumInSym, dp))) * (SIN(Angle) - SIN(AngleConj))
end do
end if
end do
end do
end do
do j = 1, NoOrbs
Norm = 0.0_dp
do i = 1, NoOrbs
Norm = Norm + (CoeffT1(i, j)**2)
end do
if (near_zero(Norm)) then
CoeffT1(j, j) = 1.0_dp
end if
end do
do j = 1, NoOrbs
do i = 1, NoOrbs
write(stdout, "(G13.5)", advance='no') CoeffT1(j, i)
end do
write(stdout, *) ""
end do
!Check normalization.
do j = 1, NoOrbs
Norm = 0.0_dp
do i = 1, NoOrbs
Norm = Norm + (CoeffT1(i, j)**2)
end do
if (abs(Norm - 1.0_dp) > 1.0e-7_dp) then
call Stop_All("EquateDiagFock", "Rotation Coefficients not normalized")
end if
end do
! Check orthogonality.
do j = 1, NoOrbs
do i = 1, NoOrbs
if (i == j) cycle
Norm = 0.0_dp
do k = 1, NoOrbs
Norm = Norm + (CoeffT1(k, j) * CoeffT1(k, i))
end do
if (abs(Norm) > 1.0e-7_dp) then
write(stdout, *) "COLUMNS: ", j, i
call Stop_All("EquateDiagFock", "RotationCoefficients not orthogonal")
end if
end do
end do
end subroutine EquateDiagFock