subroutine FindTheForce()
integer :: m, z, i, j, k, l, a, Symm, w, x, y, SymMin
real(dp) :: OffDiagForcemz, DiagForcemz, OneElForcemz, LambdaTerm1, LambdaTerm2
real(dp) :: NonDerivTerm, DerivPot
logical :: leqm, jeqm, keqm
! Running over m and z, covers all matrix elements of the force matrix (derivative
! of equation we are minimising, with respect to each translation coefficient) filling
! them in as it goes.
call set_timer(FindtheForce_time, 30)
DerivCoeff(:, :) = 0.0_dp
Force = 0.0_dp
ForceInts = 0.0_dp
OrthoForce = 0.0_dp
OffDiagForceMZ = 0
SymMin = 0
i = 0
! If the orbitals are being separated, do this whole loop twice, once for occupied and once for virtual
! i.e w = 1, 2. Otherwise do them all at once.
do w = MinOccVirt, MaxOccVirt
if (w == 1) then
SymMin = 1
MinMZ = 1
if (tSeparateOccVirt) then
MaxMZ = NoOcc
else
MaxMZ = NoOrbs
end if
else
SymMin = 9
MinMZ = NoOcc + 1
MaxMZ = NoOrbs
end if
! If we are localising the occupied and virtual orbitals separately, the above block ensures that we loop over
! first the occupied then the virtual. If we are not separating the orbitals we just run over all orbitals.
do m = MinMZ, MaxMZ
if (tStoreSpinOrbs) then
SymM = int(G1(SymLabelList2_rot(m))%sym%S)
else
SymM = int(G1(SymLabelList2_rot(m) * 2)%sym%S)
end if
do z = SymLabelCounts2_rot(1, SymM + SymMin), &
(SymLabelCounts2_rot(1, SymM + SymMin) + &
SymLabelCounts2_rot(2, SymM + SymMin) - 1)
! Find the force on a coefficient c(m,z).
OffDiagForcemz = 0.0_dp
! OffDiagForce is from any of the OffDiagMin/Max (Sqrd or not), or the double/single excitation
! max/min, as only one of these terms may be used at once.
DiagForcemz = 0.0_dp
! DiagForce includes ER localisation, and the coulomb terms <ij|ij>.
OneElForcemz = 0.0_dp
! OneElForce includes that from the one electron integrals <i|h|j> and the one particle orbital
! energies.
! DIAG TERMS
! Maximise <ii|ii>, self interaction terms.
if (tERLocalization .and. (.not. tStoreSpinOrbs)) then
DiagForcemz = DiagForcemz + (2 * ThreeIndInts01ER(z, m)) + (2 * ThreeIndInts02ER(z, m))
! Derivative of <ii|ii> only non-zero when i = m.
! each of the four terms then correspond to zeta = a, b, g, then d in the unrotated basis.
else if (tERLocalization) then
! Looking at <ij|ij> terms where j = i+1 (i.e. i is alpha of spin orbital and j is beta - or vice versa).
if (tStoreSpinOrbs) then
if (MOD(m, 2) == 0) then ! m = j
i = m - 1
DiagForcemz = DiagForcemz + ThreeIndInts01(m, i, i, z) + ThreeIndInts01(z, i, i, m) + &
ThreeIndInts01(i, m, z, i) + ThreeIndInts01(i, z, m, i)
else
j = m + 1
DiagForcemz = DiagForcemz + ThreeIndInts01(m, j, j, z) + ThreeIndInts01(z, j, j, m) + &
ThreeIndInts01(j, m, z, j) + ThreeIndInts01(j, z, m, j)
end if
else
DiagForcemz = DiagForcemz + ThreeIndInts01(m, m, m, z) + ThreeIndInts02(m, m, m, z) + &
ThreeIndInts03(m, m, m, z) + ThreeIndInts04(m, m, m, z)
end if
! First term when m = i and z = a, second when m = i and z = g.
end if
! Maximise <ij|ij>, coulomb terms, where i<j, i occ or virt, j virt only.
if (tVirtCoulombMax) then
do i = 1, NoOrbs
if (i == m) then
do j = NoOcc + 1, NoOrbs
if (j <= i) cycle ! i<j.
DiagForcemz = DiagForcemz + ThreeIndInts01(m, j, j, z) + ThreeIndInts03(m, j, j, z)
! First term for when m = i and z = a, second when m = i and z = g.
end do
end if
if ((m > NoOcc) .and. (m > i)) DiagForcemz = &
DiagForcemz + ThreeIndInts02(i, i, m, z) + ThreeIndInts04(i, i, m, z)
! This only contributes when j = m (no point in running over all j.
! First term when m = j and z = b, second when m = j and z = d.
end do
end if
! ONE ELECTRON TERMS
! Minimise |<i|h|j>|^2 where either one or bot of i and j are virtual, but i<j.
if (tHijSqrdMin) then
do j = NoOcc + 1, NoOrbs
if (m /= j) OneElForcemz = OneElForcemz + (2 * TMAT2DRot(m, j) * TMAT2DPartRot02(z, j))
! m = i and z = a.
end do
do i = NoOcc + 1, NoOrbs
if (m /= i) OneElForcemz = OneElForcemz + (2 * TMAT2DRot(i, m) * TMAT2DPartRot01(i, z))
! m = j and z = b
end do
end if
! OnePartOrbEnMax ; Maximisie sum_i [E_i - E_min]^Alpha
! where E_i = <i|h|i> + sum_j <ij||ij> and E_min is either E_LUMO (rotating virtual only) or the chemical
! potential (midway between LUMO and HOMO, when rotating all), Alpha specified in input.
! The derivative of the one part orb energies is then Alpha * NonDerivTerm * DerivPot^(Alpha-1)
! OneElIntMax ; Maximise <i|h|i>
if (tOnePartOrbEnMax .or. tOneElIntMax) then
do i = NoOcc + 1, NoOrbs
DerivPot = 0.0_dp
DerivPot = DerivPot + TMAT2DPartRot02(z, m) + TMAT2DPartRot01(m, z)
! First term when m = i and z = a, second when m = i and z = b.
! This is all that is needed for OneElIntMax
if (tOnePartOrbEnMax) then
NonDerivTerm = 0.0_dp
if (.not. (OrbEnMaxAlpha.isclose.1.0_dp)) then
! The non-derived term in the chain rule, <i|h|i> + sum_j <ij||ij> - E_min.
NonDerivTerm = NonDerivTerm + TMAT2DRot(i, i) - EpsilonMin
do j = 1, NoOcc
NonDerivTerm = NonDerivTerm + (2 * FourIndInts(i, j, i, j)) - FourIndInts(i, j, j, i)
end do
NonDerivTerm = OrbEnMaxAlpha * (NonDerivTerm**(OrbEnMaxAlpha - 1))
else
! If Alpha = 1, the NonDerivTerm will be raised to the power of 0, thus always 1.
NonDerivTerm = 1.0
end if
if (i == m) then
do j = 1, NoOcc
DerivPot = DerivPot + (2 * ThreeIndInts01(m, j, j, z)) - ThreeIndInts01(j, j, m, z) + &
(2 * ThreeIndInts03(m, j, j, z)) - ThreeIndInts03(m, j, z, j)
! First part is for when m = i and z = a, the second is for when m = i and z = g
end do
end if
! When m = j, for a particular i.
! m and z run only over virtual, and j is over occupied. m will never = j.
else
NonDerivTerm = 1.0_dp
end if
OneElForcemz = OneElForcemz + (NonDerivTerm * DerivPot)
end do
end if
! OFFDIAGTERMS
! Maximises the square of the single and double excitation integrals connected to the HF.
! I.e maximises <ij|kl> where i, j are occupied and k, l are virtual (doubles), except k may be occuppied if
! equal to i (<ij|il> singles).
! Currently this is only used for rotating virtual only, so m can only equal k or l.
if (tHFSingDoubExcMax) then
do i = 1, NoOcc
do j = 1, NoOcc
do k = NoOcc + 1, NoOrbs
if (k == m) then
do l = NoOcc + 1, NoOrbs
OffDiagForcemz = OffDiagForcemz + (2 * FourIndInts(i, j, m, l) * ThreeIndInts03(i, j, l, z))
! m = k and z = g.
end do
end if
OffDiagForcemz = OffDiagForcemz + (2 * FourIndInts(i, j, k, m) * ThreeIndInts04(i, k, j, z))
! m = l and z = d.
! Sing excitations <ij|il> where i and j are occ, l virt.
OffDiagForcemz = OffDiagForcemz + (2 * FourIndInts(i, j, i, m) * ThreeIndInts04(i, i, j, z))
! m = l
end do
end do
end do
end if
! OffDiag Sqrd/notSqrd Min/Max treats the elements <ij|kl>
! i<k and j<l.
if (tOffDiagSqrdMin .or. tOffDiagSqrdMax .or. tOffDiagMin .or. tOffDiagMax .or. tDoubExcMin) then
do l = 1, NoOrbs
if (l == m) then
leqm = .true.
else
leqm = .false.
end if
do j = 1, l - 1
if (j == m) then
jeqm = .true.
else
jeqm = .false.
end if
do k = 1, j - 1
if (k == l) cycle
if (k == m) then
keqm = .true.
else
keqm = .false.
end if
! only i with symmetry equal to j x k x l will have integrals with overall
! symmetry A1 and therefore be non-zero.
! Running across i, ThreeIndInts01 only contributes
if ((m <= k - 1) .and. (m /= j) .and. ((i /= k) .or. (j /= l))) then
if (tOffDiagSqrdMin .or. tOffDiagSqrdMax) OffDiagForcemz = OffDiagForcemz + 2 * &
(FourIndInts02(j, k, l, m) * ThreeIndInts01(k, j, l, z))
if (tOffDiagMin .or. tOffDiagMax) OffDiagForcemz = OffDiagForcemz + ThreeIndInts01(k, j, l, z)
if (tDoubExcMin) OffDiagForcemz = OffDiagForcemz + ThreeIndInts01(k, j, l, z) - &
ThreeIndInts01(l, j, k, z)
end if
if (jeqm) then
do i = 1, k - 1
if ((i /= j) .and. ((i /= k) .or. (j /= l))) then
if (tOffDiagSqrdMin .or. tOffDiagSqrdMax) OffDiagForcemz = OffDiagForcemz + 2 * &
(FourIndInts(i, j, k, l) * ThreeIndInts02(i, k, l, z))
if (tOffDiagMin .or. tOffDiagMax) &
OffDiagForcemz = OffDiagForcemz + ThreeIndInts02(i, k, l, z)
if (tDoubExcMin) OffDiagForcemz = OffDiagForcemz + (ThreeIndInts02(i, k, l, z)) - &
ThreeIndInts02(i, l, k, z)
end if
end do
end if
if (keqm) then
do i = 1, k - 1
if ((i /= j) .and. ((i /= k) .or. (j /= l))) then
if (tOffDiagSqrdMin .or. tOffDiagSqrdMax) OffDiagForcemz = OffDiagForcemz + 2 * &
(FourIndInts(i, j, k, l) * ThreeIndInts03(i, j, l, z))
if (tOffDiagMin .or. tOffDiagSqrdMax) OffDiagForcemz = OffDiagForcemz + &
ThreeIndInts03(i, j, l, z)
if (tDoubExcMin) OffDiagForcemz = OffDiagForcemz + (ThreeIndInts03(i, j, l, z)) - &
ThreeIndInts03(i, j, z, l)
end if
end do
end if
if (leqm) then
do i = 1, k - 1
if ((i /= j) .and. ((i /= k) .or. (j /= l))) then
if (tOffDiagSqrdMin .or. tOffDiagSqrdMin) OffDiagForcemz = OffDiagForcemz + 2 * &
(FourIndInts(i, j, k, l) * ThreeIndInts04(i, k, j, z))
if (tOffDiagMin .or. tOffDiagMax) &
OffDiagForcemz = OffDiagForcemz + ThreeIndInts04(i, k, j, z)
if (tDoubExcMin) OffDiagForcemz = OffDiagForcemz + (ThreeIndInts04(i, k, j, z)) - &
ThreeIndInts04(i, z, j, k)
end if
end do
end if
end do
end do
end do
end if
DerivCoeff(z, m) = (MaxMinFac * OffDiagWeight * OffDiagForcemz) + (DiagMaxMinFac * DiagWeight * DiagForcemz) + &
(OneElMaxMinFac * OneElWeight * OneElForcemz)
Force = Force + ABS(DerivCoeff(z, m))
end do
end do
end do
Force = Force / real(NoOrbs**2, dp)
! Calculate the derivatives of orthogonalisation condition.
! Have taken this out of the m and z loop to make the shake faster, but can put it back in if start using it a lot.
if (tLagrange) then
do x = MinMZ, MaxMZ
m = SymLabelList2_rot(x)
! Symmetry requirement that z must be from the same irrep as m
SymM = int(G1(m * 2)%sym%S)
do y = SymLabelCounts2_rot(1, SymM + SymMin), &
(SymLabelCounts2_rot(1, SymM + SymMin) + &
SymLabelCounts2_rot(2, SymM + SymMin) - 1)
z = SymLabelList2_rot(y)
LambdaTerm1 = 0.0_dp
LambdaTerm2 = 0.0_dp
do j = 1, NoOrbs
LambdaTerm1 = LambdaTerm1 + (Lambdas(m, j) * CoeffT1(z, j))
LambdaTerm2 = LambdaTerm2 + (Lambdas(j, m) * CoeffT1(z, j))
end do
! DerivCoeff is 'the force'. I.e. the derivative of |<ij|kl>|^2 with
! respect to each transformation coefficient. It is the values of this matrix that will tend to 0 as
! we minimise the sum of the |<ij|kl>|^2 values.
! With the Lagrange keyword this includes orthonormality conditions, otherwise it is simply the unconstrained force.
DerivCoeff(z, m) = (2 * OffDiagForcemz) - LambdaTerm1 - LambdaTerm2
OrthoForce = OrthoForce - LambdaTerm1 - LambdaTerm2
end do
end do
! If doing a Lagrange calc we also need to find the force on the lambdas to ensure orthonormality...
OrthoForce = OrthoForce / real(NoOrbs**2, dp)
DerivLambda(:, :) = 0.0_dp
do i = 1, NoOrbs
do j = 1, i
do a = 1, NoOrbs
DerivLambda(i, j) = DerivLambda(i, j) + CoeffT1(a, i) * CoeffT1(a, j)
end do
DerivLambda(j, i) = DerivLambda(i, j)
end do
end do
do i = 1, NoOrbs
DerivLambda(i, i) = DerivLambda(i, i) - 1.0_dp
end do
end if
call halt_timer(FindtheForce_Time)
end subroutine FindTheForce