subroutine Transform2ElInts()
integer :: i, j, k, l, a, b, g, d
real(dp) :: t, Temp4indints(NoRotOrbs, NoOrbs)
real(dp) :: Temp4indints02(NoRotOrbs, NoRotOrbs)
call set_timer(Transform2ElInts_time, 30)
!Zero arrays from previous transform
TwoIndInts01(:, :, :, :) = 0.0_dp
FourIndInts(:, :, :, :) = 0.0_dp
if (tNotConverged) then
TwoIndInts02(:, :, :, :) = 0.0_dp
ThreeIndInts01(:, :, :, :) = 0.0_dp
ThreeIndInts02(:, :, :, :) = 0.0_dp
ThreeIndInts03(:, :, :, :) = 0.0_dp
ThreeIndInts04(:, :, :, :) = 0.0_dp
FourIndInts02(:, :, :, :) = 0.0_dp
end if
! ************
!Transform the 1 electron, 2 index integrals (<i|h|j>).
if (tNotConverged) then
TMAT2DRot(:, :) = 0.0_dp
TMAT2DPartRot01(:, :) = 0.0_dp
TMAT2DPartRot02(:, :) = 0.0_dp
call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, &
TMAT2DTemp(:, :), NoOrbs, 0.0_dp, TMAT2DPartRot01(:, :), NoOrbs)
call dgemm('T', 'T', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, &
TMAT2DTemp(:, :), NoOrbs, 0.0_dp, TMAT2DPartRot02(:, :), NoOrbs)
call dgemm('T', 'T', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, &
TMAT2DPartRot01(:, :), NoOrbs, 0.0_dp, TMAT2DRot(:, :), NoOrbs)
end if
! **************
! Calculating the two-transformed, four index integrals.
! The untransformed <alpha beta | gamma delta> integrals are found from UMAT(UMatInd(i, j, k, l)
do b = 1, NoOrbs
do d = 1, b
Temp4indints(:, :) = 0.0_dp
call dgemm('T', 'N', NoRotOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, UMatTemp01(:, :, d, b), NoOrbs, &
0.0_dp, Temp4indints(:, :), NoRotOrbs)
! Temp4indints(i,g) comes out of here, so to transform g to k, we need the transpose of this.
Temp4indints02(:, :) = 0.0_dp
call dgemm('T', 'T', NoRotOrbs, NoRotOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, Temp4indints(:, :), NoRotOrbs, &
0.0_dp, Temp4indints02(:, :), NoRotOrbs)
! Get Temp4indits02(i,k)
do i = 1, NoRotOrbs
do k = 1, i
TwoIndInts01(d, b, k, i) = Temp4indints02(k, i)
TwoIndInts01(b, d, k, i) = Temp4indints02(k, i)
TwoIndInts01(d, b, i, k) = Temp4indints02(k, i)
TwoIndInts01(b, d, i, k) = Temp4indints02(k, i)
end do
end do
end do
end do
! These calculations are unnecessary when this routine is calculated to finalize the new orbs.
if (tNotConverged) then
do g = 1, NoOrbs
do a = 1, g
Temp4indints(:, :) = 0.0_dp
call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, UMatTemp02(:, :, a, g), NoOrbs, &
0.0_dp, Temp4indints(:, :), NoOrbs)
Temp4indints02(:, :) = 0.0_dp
call dgemm('T', 'T', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, Temp4indints(:, :), NoOrbs, &
0.0_dp, Temp4indints02(:, :), NoOrbs)
do l = 1, NoOrbs
do j = 1, l
TwoIndInts02(g, a, j, l) = Temp4indints02(j, l)
TwoIndInts02(a, g, j, l) = Temp4indints02(j, l)
TwoIndInts02(g, a, l, j) = Temp4indints02(j, l)
TwoIndInts02(a, g, l, j) = Temp4indints02(j, l)
end do
end do
end do
end do
end if
! Calculating the 3 transformed, 4 index integrals. 01 = a untransformed, 02 = b, 03 = g, 04 = d
do i = 1, NoRotOrbs
do k = 1, i
Temp4indints(:, :) = 0.0_dp
call dgemm('T', 'N', NoRotOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, TwoIndInts01(:, :, k, i), NoOrbs, &
0.0_dp, Temp4indints(:, :), NoRotOrbs)
if (tNotConverged) then
do b = 1, NoOrbs
do l = 1, NoOrbs
ThreeIndInts02(i, k, l, b) = Temp4indints(l, b)
ThreeIndInts02(k, i, l, b) = Temp4indints(l, b)
end do
end do
Temp4indints02(:, :) = 0.0_dp
call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, TwoIndInts01(:, :, k, i), NoOrbs, &
0.0_dp, Temp4indints02(:, :), NoRotOrbs)
do d = 1, NoOrbs
do j = 1, NoOrbs
ThreeIndInts04(k, i, j, d) = Temp4indints02(j, d)
ThreeIndInts04(i, k, j, d) = Temp4indints02(j, d)
end do
end do
end if
Temp4indints02(:, :) = 0.0_dp
call dgemm('T', 'T', NoRotOrbs, NoRotOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, Temp4indints(:, :), NoRotOrbs, &
0.0_dp, Temp4indints02(:, :), NoRotOrbs)
do l = 1, NoRotOrbs
do j = 1, l
FourIndInts(i, j, k, l) = Temp4indints02(j, l)
FourIndInts(i, l, k, j) = Temp4indints02(j, l)
FourIndInts(k, j, i, l) = Temp4indints02(j, l)
FourIndInts(k, l, i, j) = Temp4indints02(j, l)
if (tNotConverged) then
FourIndInts02(j, k, l, i) = Temp4indints02(j, l)
FourIndInts02(j, i, l, k) = Temp4indints02(j, l)
FourIndInts02(l, k, j, i) = Temp4indints02(j, l)
FourIndInts02(l, i, j, k) = Temp4indints02(j, l)
end if
end do
end do
end do
end do
if (tNotConverged) then
do l = 1, NoOrbs
do j = 1, l
Temp4indints(:, :) = 0.0_dp
call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, TwoIndInts02(:, :, j, l), NoOrbs, &
0.0_dp, Temp4indints(:, :), NoOrbs)
do a = 1, NoOrbs
do k = 1, NoOrbs
ThreeIndInts01(k, j, l, a) = Temp4indints(k, a)
ThreeIndInts01(k, l, j, a) = Temp4indints(k, a)
end do
end do
Temp4indints(:, :) = 0.0_dp
call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, TwoIndInts02(:, :, j, l), NoOrbs, &
0.0_dp, Temp4indints(:, :), NoOrbs)
do g = 1, NoOrbs
do i = 1, NoOrbs
ThreeIndInts03(i, l, j, g) = Temp4indints(i, g)
ThreeIndInts03(i, j, l, g) = Temp4indints(i, g)
end do
end do
end do
end do
end if
! ***************************
! Calc the potential energies for this iteration (with these transformed integrals).
! This can be sped up by merging the calculations of the potentials with the transformations, but while
! we are playing around with different potentials, it is simpler to keep these separate.
if ((.not. tReadInCoeff) .and. (.not. tUseMP2VarDenMat) .and. (.not. tFindCINatOrbs) .and. (.not. tUseHFOrbs)) then
PotEnergy = 0.0_dp
TwoEInts = 0.0_dp
PEInts = 0.0_dp
call CalcPotentials()
if (tPrintInts) call PrintIntegrals()
if ((Iteration == 0) .or. ((.not. tNotConverged) .and. (Iteration > 1))) call WriteDoubHisttofile()
if (tROHistSingExc .and. (Iteration == 0)) call WriteSingHisttofile()
! If doing Lagrange orthormalisations, find the change of the potential energy due to the orthonormality
! of the orbitals...
if (tLagrange) then
PEOrtho = 0.0_dp
do i = 1, NoOrbs
do j = 1, NoOrbs
t = 0.0_dp
do a = 1, NoOrbs
t = CoeffT1(a, i) * CoeffT1(a, j)
end do
if (i == j) t = t - 1.0_dp
PEOrtho = PEOrtho - Lambdas(i, j) * t
PotEnergy = PotEnergy - Lambdas(i, j) * t
end do
end do
end if
end if
call halt_timer(Transform2ElInts_Time)
end subroutine Transform2ElInts