CalcPotentials Subroutine

    subroutine CalcPotentials()

        ! Only temporarily like this, can tidy it up majorly

        integer :: i, j, k, l, Starti, Finishi
        real(dp) :: MaxTerm

        l = 0
        if (tERLocalization .and. (.not. tStoreSpinOrbs)) then
            ERPotEnergy = 0.0_dp
            if (tRotateVirtOnly) then
                Starti = NoOcc + 1
                Finishi = NoOrbs
            else if (tRotateOccOnly) then
                Starti = 1
                Finishi = NoOcc
            else
                Starti = 1
                Finishi = NoOrbs
            end if
            CoulPotEnergy = 0.0_dp
            OffDiagPotEnergy = 0.0_dp
            do i = Starti, Finishi
                ERPotEnergy = ERPotEnergy + FourIndIntsER(i)
                if (FourIndIntsER(i) < 0) then
                    call neci_flush(stdout)
                    call Stop_All('CalcPotentials', 'A <ii|ii> value is less than 0.')
                end if
                PotEnergy = PotEnergy + FourIndIntsER(i)
                TwoEInts = TwoEInts + FourIndIntsER(i)
                PEInts = PEInts + FourIndIntsER(i)
            end do
        else if (tERLocalization) then
            ERPotEnergy = 0.0_dp
            PotEnergy = 0.0_dp
            if (tRotateVirtOnly) then
                Starti = NoOcc + 1
                Finishi = NoOrbs
            else if (tRotateOccOnly) then
                Starti = 1
                Finishi = NoOcc
            else
                Starti = 1
                Finishi = NoOrbs
            end if
            do i = Starti, Finishi
                if (tStoreSpinOrbs) then
                    if (MOD(i, 2) == 0) then
                        j = i - 1
                    else
                        j = i + 1
                    end if
                    ERPotEnergy = ERPotEnergy + FourIndInts(i, j, i, j)
                    if ((FourIndInts(i, j, i, j) < 0) .or. (FourIndInts(j, i, j, i) < 0)) then
                        call neci_flush(stdout)
                        call Stop_All('CalcPotentials', 'A <ii|ii> value is less than 0.')
                    end if
                    PotEnergy = PotEnergy + FourIndInts(i, j, i, j)
                    TwoEInts = TwoEInts + FourIndInts(i, j, i, j)
                    PEInts = PEInts + FourIndInts(i, j, i, j)
                else
                    ERPotEnergy = ERPotEnergy + FourIndInts(i, i, i, i)
                    if ((FourIndInts(i, i, i, i) < 0)) then
                        call neci_flush(stdout)
                        call Stop_All('CalcPotentials', 'A <ii|ii> value is less than 0.')
                    end if
                    PotEnergy = PotEnergy + FourIndInts(i, i, i, i)
                    TwoEInts = TwoEInts + FourIndInts(i, i, i, i)
                    PEInts = PEInts + FourIndInts(i, i, i, i)
                end if
            end do
        end if

        if (tOffDiagSqrdMin .or. tOffDiagSqrdMax .or. tOffDiagMin .or. tOffdiagMax) then
            do l = 1, NoOrbs
                do j = 1, l - 1
                    do k = 1, j - 1
                        do i = 1, k - 1
                            if (tOffDiagSqrdMin .or. tOffDiagSqrdMax) then
                                if (((i /= j) .and. (j /= l)) .and. ((i /= k) .or. (j /= l))) then
                                    PotEnergy = PotEnergy + (FourIndInts(i, j, k, l)**2)
                                    TwoEInts = TwoEInts + (FourIndInts(i, j, k, l)**2)
                                    PEInts = PEInts + (FourIndInts(i, j, k, l)**2)
                                end if
                            end if
                            if (tOffDiagMin .or. tOffDiagMax) then
                                if (.not. ((k == i) .or. (j == l))) then
                                    PotEnergy = PotEnergy + FourIndInts(i, j, k, l)
                                    TwoEInts = TwoEInts + FourIndInts(i, j, k, l)
                                    PEInts = PEInts + FourIndInts(i, j, k, l)
                                end if
                            end if
                        end do
                    end do
                end do
            end do
        end if

        if (tDoubExcMin) then
            do i = 1, NoOrbs
                do j = 1, NoOrbs
                    do k = 1, i - 1
                        if ((k == l) .and. (k == i)) cycle
                        do l = 1, j - 1
                            if ((j == k) .and. (j == l)) cycle
                            if ((j == k) .and. (j == i)) cycle
                            if ((j == l) .and. (j == i)) cycle
                            PotEnergy = PotEnergy + (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
                            TwoEInts = TwoEInts + (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
                            PEInts = PEInts + (FourIndInts(i, j, k, l)) - FourIndInts(i, j, l, k)
                        end do
                    end do
                end do
            end do
        end if

        if (tOnePartOrbEnMax .or. tOneElIntMax) then
            do i = NoOcc + 1, NoOrbs
                MaxTerm = 0.0_dp
                MaxTerm = TMAT2DRot(i, i)
                if (tOnePartOrbEnMax) then
                    do j = 1, NoOcc
                        MaxTerm = MaxTerm + (2 * FourIndInts(i, j, i, j)) - FourIndInts(i, j, j, i)
                    end do
                    MaxTerm = MaxTerm - EpsilonMin
                    MaxTerm = MaxTerm**OrbEnMaxAlpha
                end if
                PotEnergy = PotEnergy + MaxTerm
            end do
        end if

        if (tHijSqrdMin) then
            HijSqrdPotEnergy = 0.0_dp
            do i = NoOcc + 1, NoOrbs
                do j = NoOcc + 1, NoOrbs
                    if (j > i) then
                        PotEnergy = PotEnergy + (TMAT2DRot(i, j)**2)
                        HijSqrdPotEnergy = HijSqrdPotEnergy + (TMAT2DRot(i, j)**2)
                    end if
                end do
            end do
        end if

        if (tVirtCoulombMax) then
            ERPotEnergy = 0.0_dp
            ijOccVirtPotEnergy = 0.0_dp
            do i = 1, NoOrbs
                if (i <= NoOcc) then
                    do j = NoOcc + 1, NoOrbs
                        ijOccVirtPotEnergy = ijOccVirtPotEnergy + FourIndInts(i, j, i, j)
                    end do
                end if
                if (i > NoOcc) then
                    ERPotEnergy = ERPotEnergy + FourIndInts(i, i, i, i)
                    do j = NoOcc + 1, NoOrbs
                        if (j <= i) cycle
                        PotEnergy = PotEnergy + FourIndInts(i, j, i, j)
                        TwoEInts = TwoEInts + FourIndInts(i, j, i, j)
                    end do
                end if
            end do
        end if

        if (tHFSingDoubExcMax) then
            do i = 1, NoOcc
                do j = 1, NoOcc
                    do k = NoOcc + 1, NoOrbs
                        do l = NoOcc + 1, NoOrbs
                            PotEnergy = PotEnergy + (FourIndInts(i, j, k, l)**2)
                        end do

                        ! Sing excitations <ij|ik> where i and j are occ, k virt.
                        PotEnergy = PotEnergy + (FourIndInts(i, j, i, k)**2)
                    end do
                end do
            end do
        end if

    end subroutine CalcPotentials