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