subroutine BrokenSymNO(evalues, occ_numb_diff)
! This rouine finds natural orbitals (NOs) whose occupation
! numbers differ by a small relative threshold (occ_numb_diff) and
! rotates them by calling the Rotate2Orbs routine in order to
! break symmetry and maximally localise the NOs
! We'd like to keep both the original NOs (Natural Orbitals)
! and the broken maximally localised NOs for checking.
! This is not very (time-) efficient at the moment.
use IntegralsData, only: umat
use LoggingData, only: tBreakSymNOs, local_cutoff, tagRotNOs, RotNOs
use LoggingData, only: rottwo, rotthree, rotfour
use MemoryManager, only: LogMemDealloc
use Parallel_neci, only: iProcIndex
use rdm_data, only: tRotatedNOs
use SystemData, only: nel, tStoreSpinOrbs
use UMatCache, only: UMatInd
real(dp), intent(in) :: evalues(:)
real(dp), intent(in) :: occ_numb_diff
real(dp) :: diffnorm, SumDiag, sum_old, sum_new, selfint_old
real(dp), allocatable :: one_rdm(:, :), trans_2orbs_coeffs(:, :)
real(dp), allocatable :: selfint(:)
integer, allocatable :: rotate_list(:, :), rotorbs(:, :)
integer, allocatable :: sym_list(:)
integer :: l1, l2, l3, l4, l5, l6, m, n
integer :: iumat, jumat
logical :: partnerfound, localdelocal
allocate(one_rdm(NoOrbs, NoOrbs))
allocate(sym_list(NoOrbs))
allocate(trans_2orbs_coeffs(2, 2))
allocate(rotorbs(6, 2))
allocate(rotate_list((NoOrbs * (NoOrbs - 1)), 4))
allocate(selfint(NoOrbs))
if (iProcIndex == 0) then
! No symmetry ordering is applied.
do l1 = 1, NoOrbs
sym_list(l1) = l1
end do
! Normalisation.
SumDiag = sum(evalues)
if (tOpenShell) then
diffnorm = SumDiag / dble(NEl)
else
diffnorm = 2.0_dp * (SumDiag / dble(NEl))
end if
trotatedNOs = .true.
if (tStoreSpinorbs) then
call Stop_all("BrokenSymNO", "Broken symmetry NOs currently not implemented for UHF")
end if
write(stdout, *) '------------------------------------------------------------------------------'
write(stdout, *) 'Localising NOs whose occupation numbers differ by less than threshold'
write(stdout, *) '------------------------------------------------------------------------------'
if (tBreakSymNOs) then
write(stdout, *) 'Rotating specified NOs'
else
write(stdout, *) 'Threshold for orbitals to rotate:', occ_numb_diff
end if
! Self-interactions.
selfint(:) = 0.0_dp
do l1 = 1, NoOrbs
selfint(l1) = Umat(UmatInd(l1, l1, l1, l1))
end do
write(stdout, *) 'Self-interactions for NOs:'
do l1 = 1, NoOrbs
write(stdout, '(I3,3X,G25.12)') l1, selfint(l1)
end do
write(stdout, *) 'Sum of NO selfinteractions:', sum(selfint)
selfint_old = sum(selfint)
! If the NOs to be rotated are specified in the input file.
if (tBreakSymNOs) then
rotate_list(:, :) = 0
rotorbs(:, :) = 0
m = 0
do l1 = 2, (2 * rottwo), 2
m = m + 1
rotate_list(m, 1) = RotNOs(l1 - 1)
rotate_list(m, 2) = RotNOs(l1)
end do
do l1 = ((2 * rottwo) + 3), ((2 * rottwo) + (3 * rotthree)), 3
m = m + 1
rotate_list(m, 1) = RotNOs(l1 - 2)
rotate_list(m, 2) = RotNOs(l1 - 1)
rotate_list(m, 3) = RotNOs(l1)
end do
do l1 = ((2 * rottwo) + (3 * rotthree) + 4), ((2 * rottwo) + (3 * rotthree) + (4 * rotfour)), 4
m = m + 1
rotate_list(m, 1) = RotNOs(l1 - 3)
rotate_list(m, 2) = RotNOs(l1 - 2)
rotate_list(m, 3) = RotNOs(l1 - 1)
rotate_list(m, 4) = RotNOs(l1)
end do
else
! If the threshold is used to generate a list of NOs to be
! rotated.
! Generate the list of orbitals which are rotated amongst each
! other.
rotate_list(:, :) = 0
rotorbs(:, :) = 0
! Need to account for spatial and spin orbital representations
! since orbitals of different spin cannot be mixed.
! List contains the NOs which are rotated.
! It can deal with a maximum of four NOs which are mixed.
m = 0
n = 1
do l1 = 1, NoOrbs
if ((m /= 0) .and. (l1 <= rotate_list(m, n))) cycle
partnerfound = .false.
n = 1
do l2 = (l1 + 1), NoOrbs
if ((abs((evalues(l1) / diffnorm) - (evalues(l2) / diffnorm)) / abs((evalues(l2) / diffnorm)))&
& < occ_numb_diff) then
if (.not. partnerfound) then
m = m + 1
n = n + 1
rotate_list(m, 1) = l1
rotate_list(m, 2) = l2
partnerfound = .true.
else if (partnerfound) then
n = n + 1
! this is for up to 2-fold degenearcy
! if (n.gt.2) then
! n = 2
! write(stdout,*) '***Warning***'
! write(stdout,*) 'Threshold generated more than 2-fold degeneracy'
! write(stdout,*) 'NOs around:',l2
! cycle ! don't want to rotate more than 2 orbitals
! end if
! this is for up to four-fold degeneracy
if (n > 4) then
n = 4
write(stdout, *) '***Warning***'
write(stdout, *) 'Threshold generated more than 4-fold degeneracy'
write(stdout, *) 'NOs around:', l2
cycle ! don't want to rotate more than 4 orbitals
end if
rotate_list(m, n) = l2
end if
end if
end do
end do
end if
write(stdout, *) 'The following pairs of orbitals will be rotated:'
do l1 = 1, m
write(stdout, '(I3,3X,4(I3))') l1, rotate_list(l1, :)
end do
one_rdm(:, :) = 0.0_dp
do l1 = 1, NoOrbs
one_rdm(l1, l1) = 1.0_dp
end do
! Rotate two-fold degenerate pairs first.
do l1 = 1, m
! If only two orbitals have the same occupation numbers.
if (rotate_list(l1, 3) == 0) then
write(stdout, '(A20,4(I3))') 'Rotating NOs:', rotate_list(l1, :)
iumat = rotate_list(l1, 1)
jumat = rotate_list(l1, 2)
if (jumat <= local_cutoff) then
localdelocal = .false.
else if (jumat > local_cutoff) then
localdelocal = .true.
end if
call Rotate2Orbs(iumat, jumat, trans_2orbs_coeffs, selfint(iumat),&
&selfint(jumat), localdelocal)
! The new NOs are
! phi_{i'} = cos a p_{i} + sin a p_{j}
! phi_{j'} = -sin a p_{i} + cos a p_{j}
one_rdm(iumat, iumat) = trans_2orbs_coeffs(1, 1)
one_rdm(jumat, iumat) = trans_2orbs_coeffs(2, 1)
one_rdm(iumat, jumat) = trans_2orbs_coeffs(1, 2)
one_rdm(jumat, jumat) = trans_2orbs_coeffs(2, 2)
write(stdout, *) 'Sum of rotated NO self-interactions:', sum(selfint)
selfint_old = sum(selfint)
end if
end do
! Transform integral corresponding to rotated NO.
! These are required for not printing out RPFCIDUMP or
! BSFCIDUMP every time.
trotatedNOs = .true.
call Transform2ElIntsMemSave_RDM(one_rdm, sym_list)
call RefillUMATandTMAT2D_RDM(one_rdm, sym_list)
! If three orbitals are degenerate.
do l1 = 1, m
if ((rotate_list(l1, 3) /= 0) .and. (rotate_list(l1, 4) == 0)) then
sum_new = sum(selfint)
rotorbs(1, 1) = 1
rotorbs(1, 2) = 2
rotorbs(2, 1) = 1
rotorbs(2, 2) = 3
rotorbs(3, 1) = 2
rotorbs(3, 2) = 3
! These have to be done self-consistently since all
! three orbitals can intermix.
do
sum_old = sum_new
write(stdout, '(A20,4(I3))') 'Rotating NOs:', rotate_list(l1, :)
do l3 = 1, 3
iumat = rotate_list(l1, rotorbs(l3, 1))
jumat = rotate_list(l1, rotorbs(l3, 2))
one_rdm = 0.0_dp
do l4 = 1, NoOrbs
if ((l4 /= iumat) .and. (l4 /= jumat)) then
one_rdm = 1.0_dp
end if
end do
if (jumat <= local_cutoff) then
localdelocal = .false.
else if (jumat > local_cutoff) then
localdelocal = .true.
end if
call Rotate2Orbs(iumat, jumat,&
&trans_2orbs_coeffs, selfint(iumat), selfint(jumat)&
&, localdelocal)
! The new NOs are
! phi_{i'} = cos a p_{i} + sin a p_{j}
! phi_{j'} = -sin a p_{i} + cos a p_{j}
one_rdm(iumat, iumat) = trans_2orbs_coeffs(1, 1)
one_rdm(jumat, iumat) = trans_2orbs_coeffs(2, 1)
one_rdm(iumat, jumat) = trans_2orbs_coeffs(1, 2)
one_rdm(jumat, jumat) = trans_2orbs_coeffs(2, 2)
! Transform integral corresponding to rotated NO.
! These are required for not printing out
! RPFCIDUMP or BSFCIDUMP every time.
trotatedNOs = .true.
call Transform2ElIntsMemSave_RDM(one_rdm, sym_list)
call RefillUMATandTMAT2D_RDM(one_rdm, sym_list)
end do
! Check for convergence.
sum_new = sum(selfint)
write(stdout, '(A50,2G20.12)') 'Current and previous selfinteraction:',&
&sum_new, sum_old
if (abs(sum_new - sum_old) < 1e-12_dp) then
exit
end if
end do
write(stdout, *) 'Sum of rotated NO self-interactions:', sum(selfint)
selfint_old = sum(selfint)
else if ((rotate_list(l1, 3) /= 0) .and. (rotate_list(l1, 4) /= 0)) then
sum_new = sum(selfint)
rotorbs(1, 1) = 1
rotorbs(1, 2) = 2
rotorbs(2, 1) = 3
rotorbs(2, 2) = 4
rotorbs(3, 1) = 1
rotorbs(3, 2) = 3
rotorbs(4, 1) = 2
rotorbs(4, 2) = 4
rotorbs(5, 1) = 1
rotorbs(5, 2) = 4
rotorbs(6, 1) = 2
rotorbs(6, 2) = 3
! These have to be done self-consistently since all
! three orbitals can intermix.
do
sum_old = sum_new
write(stdout, '(A20,4(I3))') 'Rotating NOs:', rotate_list(l1, :)
do l3 = 1, 3
one_rdm(:, :) = 0.0_dp
do l4 = 1, NoOrbs
if ((l4 /= rotate_list(l1, 1)) .and. (l4 /= rotate_list(l1, 2))&
& .and. (l4 /= rotate_list(l1, 3)) .and. (l4 /= rotate_list(l1, 4))) then
one_rdm(l4, l4) = 1.0_dp
end if
end do
! Rotate these two independently.
do l5 = 0, 1
iumat = rotate_list(l1, rotorbs(((2 * l3) - l5), 1))
jumat = rotate_list(l1, rotorbs(((2 * l3) - l5), 2))
if (jumat <= local_cutoff) then
localdelocal = .false.
else if (jumat > local_cutoff) then
localdelocal = .true.
end if
call Rotate2Orbs(iumat, jumat, trans_2orbs_coeffs, &
selfint(iumat), selfint(jumat), localdelocal)
! The new NOs are
! phi_{i'} = cos a p_{i} + sin a p_{j}
! phi_{j'} = -sin a p_{i} + cos a p_{j}
one_rdm(iumat, iumat) = trans_2orbs_coeffs(1, 1)
one_rdm(jumat, iumat) = trans_2orbs_coeffs(2, 1)
one_rdm(iumat, jumat) = trans_2orbs_coeffs(1, 2)
one_rdm(jumat, jumat) = trans_2orbs_coeffs(2, 2)
end do
! Transform integral corresponding to rotated NO.
! These are required for not printing out
! RPFCIDUMP or BSFCIDUMP every time.
trotatedNOs = .true.
call Transform2ElIntsMemSave_RDM(one_rdm, sym_list)
call RefillUMATandTMAT2D_RDM(one_rdm, sym_list)
end do
! Check for convergence.
sum_new = sum(selfint)
write(stdout, "(A50,2G20.12)") 'Current and previous selfinteractions:',&
&sum_new, sum_old
if (abs(sum_new - sum_old) < 1e-12_dp) then
exit
end if
end do
write(stdout, *) 'Sum of rotated NO self-interactions:', sum(selfint)
selfint_old = sum(selfint)
end if
end do
write(stdout, *) 'Final self-interactions for rotated NOs:'
do l1 = 1, NoOrbs
write(stdout, '(I3,3X,G25.12)') l1, selfint(l1)
end do
write(stdout, *) 'Sum of rotated NO self-interactions:', sum(selfint)
write(stdout, *) '------------------------------------------------------'
write(stdout, *) 'Writing out BSFCIDUMP...'
write(stdout, *) '------------------------------------------------------'
call PrintROFCIDUMP_RDM("BSFCIDUMP")
end if
deallocate(one_rdm)
deallocate(sym_list)
deallocate(trans_2orbs_coeffs)
deallocate(rotate_list)
deallocate(rotorbs)
deallocate(selfint)
if (tBreakSymNOs) then
deallocate(RotNOs)
call LogMemDealloc('BrokenSymNO', tagRotNOs)
end if
end subroutine BrokenSymNO