subroutine Rotate2Orbs(i, j, trans_2orbs_coeffs, selfintorb1, selfintorb2, localdelocal)
! This routine takes two orbitals i,j, and rotates them in order to
! maximally localise these. It employs an Edminston-Ruedenberg type
! localisation which maximises the self-interaction
! \sum_{i=1}^{2} \sum_{r,s,u,v} (c_{ir})*(c_{is})*c_{iu}c_{iv} <p_{i}p_{i}|u|p_{i}p_{i}>
! where p_{i} are the original NOs.
! The coefficients c are given by the following matrix:
! c = cos a sin a
! -sin a cos a
! Then angle a is found by differentiating and setting it equal to 0
! which gives the following analytical expression of the form
! tan a = -x/y
! where x and y are sums of the original NO four index integrals.
use IntegralsData, only: umat
use UMatCache, only: UMatInd
real(dp), allocatable, intent(inout) :: trans_2orbs_coeffs(:, :)
real(dp), intent(inout) :: selfintorb1, selfintorb2
real(dp) :: alpha2(17)
real(dp) :: selfinteractions(17)
real(dp) :: coeffcos, coeffsin, maxint
integer :: maxangle(1)
integer :: indicesij(2)
integer, intent(in) :: i, j
integer :: l1, l2, l3, l4, l5
logical, intent(in) :: localdelocal
indicesij(1) = i
indicesij(2) = j
trans_2orbs_coeffs(:, :) = 0.0_dp
! Umat(UMatInd(i,j,k,l)) contains the four-index integrals
! <ij|kl> (physical notation) in the NO basis
coeffcos = Umat(UmatInd(i, i, i, j)) + Umat(UmatInd(i, i, j, i)) + Umat(UmatInd(i, j, i, i)) &
& - Umat(UmatInd(i, j, j, j)) + Umat(UmatInd(j, i, i, i)) - Umat(UmatInd(j, i, j, j)) &
& - Umat(UmatInd(j, j, i, j)) - Umat(UmatInd(j, j, j, i))
coeffsin = -Umat(UmatInd(i, i, i, i)) + Umat(UmatInd(i, i, j, j)) + Umat(UmatInd(i, j, i, j)) &
& + Umat(UmatInd(i, j, j, i)) + Umat(UmatInd(j, i, i, j)) + Umat(UmatInd(j, i, j, i)) &
& + Umat(UmatInd(j, j, i, i)) - Umat(UmatInd(j, j, j, j))
! atan return a value in [-pi/2,pi/2]
! because of the 4*alpha in the equation there are 8 distinct solutions
! i.e. in the range 0,2*pi
! i.e. possible solutions are separated by (2*pi/8)=pi/4
! for safety 16 solutions are evaluated.
alpha2(9) = atan((-coeffcos / coeffsin))
alpha2(9) = alpha2(9) / 4.0_dp
do l1 = 8, 1, -1
alpha2(l1) = alpha2(l1 + 1) - (pi / 4.0_dp)
end do
do l1 = 10, 17
alpha2(l1) = alpha2(l1 - 1) + (pi / 4.0_dp)
end do
! second derivatives to find maximum (necessary since the minimum, i.e. fully delocalised
! orbitals satisfy the same conditions
!secondderiv(1) = (4.0_dp*coeffsin*cos(4.0_dp*alpha(1))) - (4.0_dp*coeffcos*sin(4.0_dp*alpha(1)))
!secondderiv(2) = (4.0_dp*coeffsin*cos(4.0_dp*alpha(2))) - (4.0_dp*coeffcos*sin(4.0_dp*alpha(2)))
! Compute selfinteractions to check which one is largest.
! This is a better measure than the second derivatives.
selfinteractions(:) = 0.0_dp
do l1 = 1, 17
trans_2orbs_coeffs(1, 1) = cos(alpha2(l1))
trans_2orbs_coeffs(2, 1) = sin(alpha2(l1))
trans_2orbs_coeffs(1, 2) = -sin(alpha2(l1))
trans_2orbs_coeffs(2, 2) = cos(alpha2(l1))
do l2 = 1, 2
do l3 = 1, 2
do l4 = 1, 2
do l5 = 1, 2
selfinteractions(l1) = selfinteractions(l1) + trans_2orbs_coeffs(l2, 1)&
&* trans_2orbs_coeffs(l3, 1) *&
&trans_2orbs_coeffs(l4, 1) * trans_2orbs_coeffs(l5, 1) *&
&Umat(UmatInd(indicesij(l2), indicesij(l3), indicesij(l4), indicesij(l5)))
end do
end do
end do
end do
do l2 = 1, 2
do l3 = 1, 2
do l4 = 1, 2
do l5 = 1, 2
selfinteractions(l1) = selfinteractions(l1) + trans_2orbs_coeffs(l2, 2)&
&* trans_2orbs_coeffs(l3, 2) *&
&trans_2orbs_coeffs(l4, 2) * trans_2orbs_coeffs(l5, 2) *&
&Umat(UmatInd(indicesij(l2), indicesij(l3), indicesij(l4), indicesij(l5)))
end do
end do
end do
end do
end do
! Choose the angle which maximises the self interactions.
if (.not. localdelocal) then
! Maximally delocalised.
maxangle = minloc(selfinteractions)
maxint = minval(selfinteractions)
else if (localdelocal) then
! Maximally localised.
maxangle = maxloc(selfinteractions)
maxint = maxval(selfinteractions)
end if
! Return transformatin coefficients.
trans_2orbs_coeffs(1, 1) = cos(alpha2(maxangle(1)))
trans_2orbs_coeffs(2, 1) = sin(alpha2(maxangle(1)))
trans_2orbs_coeffs(1, 2) = -sin(alpha2(maxangle(1)))
trans_2orbs_coeffs(2, 2) = cos(alpha2(maxangle(1)))
! New self-interactions for transformed orbitals.
selfintorb1 = 0.0_dp
selfintorb2 = 0.0_dp
do l2 = 1, 2
do l3 = 1, 2
do l4 = 1, 2
do l5 = 1, 2
selfintorb1 = selfintorb1 + trans_2orbs_coeffs(l2, 1)&
&* trans_2orbs_coeffs(l3, 1) *&
&trans_2orbs_coeffs(l4, 1) * trans_2orbs_coeffs(l5, 1) *&
&Umat(UmatInd(indicesij(l2), indicesij(l3), indicesij(l4), indicesij(l5)))
selfintorb2 = selfintorb2 + trans_2orbs_coeffs(l2, 2)&
&* trans_2orbs_coeffs(l3, 2) *&
&trans_2orbs_coeffs(l4, 2) * trans_2orbs_coeffs(l5, 2) *&
&Umat(UmatInd(indicesij(l2), indicesij(l3), indicesij(l4), indicesij(l5)))
end do
end do
end do
end do
end subroutine Rotate2Orbs