Transform2ElIntsERlocal Subroutine

    subroutine Transform2ElIntsERlocal()

        integer :: i, j, a, b, g, d, m
        real(dp) :: t, Temp4indints(NoOrbs, NoOrbs)
        real(dp) :: Temp4indints02(NoOrbs)

        call set_timer(Transform2ElInts_time, 30)

        ! Zero arrays from previous transform.

        TwoIndIntsER(:, :, :) = 0.0_dp

        ThreeIndInts01ER(:, :) = 0.0_dp
        ThreeIndInts02ER(:, :) = 0.0_dp

        FourIndIntsER(:) = 0.0_dp

! **************
! Calculating the two-transformed, four index integrals.

! The untransformed <alpha beta | gamma delta> integrals are found from UMAT(UMatInd(i, j, k, l)

        do d = 1, NoOrbs
            do b = 1, d
                Temp4indints(:, :) = 0.0_dp
                Temp4indints02(:) = 0.0_dp
                call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, UMATTemp01(:, :, b, d), &
                           NoOrbs, 0.0_dp, Temp4indints(:, :), NoOrbs)
                ! a -> m. Temp4indints(m,g) comes out of here.
                ! Want to transform g to m as well.

                do m = 1, NoOrbs
                    do g = 1, NoOrbs
                        Temp4indints02(m) = Temp4indints02(m) + (Temp4indints(m, g) * CoeffT1(g, m))
                    end do
                end do
                ! Now have Temp4indints(m,m) for each b and d.

                do m = 1, NoOrbs
                    TwoIndIntsER(b, d, m) = Temp4indints02(m)
                    TwoIndIntsER(d, b, m) = Temp4indints02(m)
                end do
            end do
        end do

! Now want to transform g to get one of the 3-transformed 4-index integrals <a m | m m>.
! These can be stored in 2-D arrays, as they can be specified by only m and z.

        do m = 1, NoOrbs
            do b = 1, NoOrbs
                do d = 1, NoOrbs
                    ThreeIndInts01ER(b, m) = ThreeIndInts01ER(b, m) + (TwoIndIntsER(b, d, m) * CoeffT1(d, m))
                end do
            end do
        end do
        ! ThreeIndInts01ER(z,m) is where z is alpha (a).

        TwoIndIntsER(:, :, :) = 0.0_dp
        do d = 1, NoOrbs
            do b = 1, NoOrbs
                Temp4indints(:, :) = 0.0_dp
                Temp4indints02(:) = 0.0_dp
                call dgemm('T', 'N', NoOrbs, NoOrbs, NoOrbs, 1.0_dp, CoeffT1(:, :), NoOrbs, UMATTemp01(:, :, b, d), &
                           NoOrbs, 0.0_dp, Temp4indints(:, :), NoOrbs)
                ! a -> m. Temp4indints(m,g) comes out of here.
                ! Want to transform g to m as well.

                do m = 1, NoOrbs
                    do g = 1, NoOrbs
                        Temp4indints02(m) = Temp4indints02(m) + (Temp4indints(m, g) * CoeffT1(g, m))
                    end do
                end do
                ! Now have Temp4indints(m,m) for each a and g.

                do m = 1, NoOrbs
                    TwoIndIntsER(b, d, m) = Temp4indints02(m)
                    TwoIndIntsER(d, b, m) = Temp4indints02(m)
                end do
            end do
        end do

! Now want to transform g to get one of the 3-transformed 4-index integrals <a m | m m>.
! These can be stored in 2-D arrays, as they can be specified by only m and z.

        do m = 1, NoOrbs
            do b = 1, NoOrbs
                do d = 1, NoOrbs
                    ThreeIndInts02ER(b, m) = ThreeIndInts02ER(b, m) + (TwoIndIntsER(b, d, m) * CoeffT1(d, m))
                end do
            end do
        end do
        ! ThreeIndInts02ER(z,m) is where z is beta (b).

! Find the <ii|ii> integrals, to calculate the potential energy.

        do m = 1, NoOrbs
            do a = 1, NoOrbs
                FourIndIntsER(m) = FourIndIntsER(m) + (ThreeIndInts01ER(a, m) * CoeffT1(a, m))
            end do
        end do

! ***************************
! 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.

        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

        call halt_timer(Transform2ElInts_Time)

    end subroutine Transform2ElIntsERlocal