Transform2ElInts Subroutine

    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