I have some Python code containing a couple of for loops, which I would like to optimize by using low-level functions. My first approach is to use numpy.einsum(), but the results don't match the for-loop implementation:## first termT11 = 0.125 * numpy.einsum('abcd,efab,cdef,abcd,efcd',    g[:Nelec,:Nelec,Nelec:,Nelec:], # <ab||rs> -> abcd    g[:Nelec,:Nelec,:Nelec,:Nelec], # <cd||ab> -> efab    g[Nelec:,Nelec:,:Nelec,:Nelec], # <rs||cd> -> cdef    D,                              # e_a + e_b - e_r - e_s -> abcd    D,                              # e_c + e_d - e_r - e_s -> efcd)T12 = 0for a in range(Nelec):    for b in range(Nelec):        for c in range(Nelec):            for d in range(Nelec):                for r in range(Nvirt):                    for s in range(Nvirt):                        # nominator                        nom  = g[a,b,Nelec+r,Nelec+s]                        nom *= g[c,d,a,b]                        nom *= g[Nelec+r,Nelec+s,c,d]                        # denominator                        denom = D[a,b,r,s] * D[c,d,r,s]                        T12 += nom*denomT12 *= 0.125## second termT21 = 0.125 * numpy.einsum('abcd,cdef,efab,abcd,abef',    g[:Nelec,:Nelec,Nelec:,Nelec:], # <ab||rs> -> abcd    g[Nelec:,Nelec:,Nelec:,Nelec:], # <rs||tu> -> cdef    g[Nelec:,Nelec:,:Nelec,:Nelec], # <tu||ab> -> efab    D,                              # e_a + e_b - e_r - e_s -> abcd    D,                              # e_a + e_b - e_t - e_u -> abef)T22 = 0for a in range(Nelec):    for b in range(Nelec):        for r in range(Nvirt):            for s in range(Nvirt):                for t in range(Nvirt):                    for u in range(Nvirt):                        # nominator                        nom  = g[a,b,Nelec+r,Nelec+s]                        nom *= g[Nelec+r,Nelec+s,Nelec+t,Nelec+u]                        nom *= g[Nelec+t,Nelec+u,a,b]                        # denominator                        denom = D[a,b,r,s] * D[a,b,t,u]                        T22 += nom*denomT22 *= 0.125## third termT31 = numpy.einsum('abcd,edfb,cfae,abcd,aecf',    g[:Nelec,:Nelec,Nelec:,Nelec:], # <ab||rs> -> abcd    g[:Nelec,Nelec:,Nelec:,:Nelec], # <cs||tb> -> edfb    g[Nelec:,Nelec:,:Nelec,:Nelec], # <rt||ac> -> cfae    D,                              # e_a + e_b - e_r - e_s -> abcd    D,                              # e_a + e_c - e_r - e_t -> aecf)T32 = 0for a in range(Nelec):    for b in range(Nelec):        for c in range(Nelec):            for r in range(Nvirt):                for s in range(Nvirt):                    for t in range(Nvirt):                        # nominator                        nom  = g[a,b,Nelec+r,Nelec+s]                        nom *= g[c,Nelec+s,Nelec+t,b]                        nom *= g[Nelec+r,Nelec+t,a,c]                        # denominator                        denom = D[a,b,r,s] * D[a,c,r,t]                        T32 += nom*denomprint(T11, T12)print(T21, T22)print(T31, T32)print(T11+T21+T31, T12+T22+T32)For my test case I get the following output: 0.016038311274  0.016038311274 0.014256694984  0.014256694984-0.044483006105 -0.044483006105-0.014187999845 -0.014187999845The first column contains the einsum() results, the second column is for the correctly working for loop implementation. So my question is: Why does einsum() not give the same results (except for the second term $T_2$)? Is there some error in my einsum() implementation or is there a conceptual problem here?Edit: After fixing one error in T11 and another error in T31 both methods yield same results.Not necessarily  important, but the equations are for the Mller-Plesset perturbation theory third-order energy correction:\begin{equation} E_0^{(3)} = T_1 + T_2 + T_3 \\ T_1 = \frac{1}{8} \sum_{abcdrs} \frac{\langle ab||rs\rangle\langle cd||ab\rangle\langle rs||cd\rangle}{(\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s)(\epsilon_c+\epsilon_d-\epsilon_r-\epsilon_s)} \\ T_2 = \frac{1}{8} \sum_{abrstu} \frac{\langle ab||rs\rangle\langle rs||tu\rangle\langle tu||ab\rangle}{(\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s)(\epsilon_c+\epsilon_d-\epsilon_t-\epsilon_u)} \\ T_3 = \sum_{abcrst} \frac{\langle ab||rs\rangle\langle cs||tb\rangle\langle rt||ac\rangle}{(\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s)(\epsilon_a+\epsilon_c-\epsilon_r-\epsilon_t)}\end{equation}In the code I use:g[i,j,k,l] = $\langle ij||kl\rangle$ being the anti-symmetrized two electron repulsion integralsD[a,b,r,s] = $\frac{1}{\epsilon_i+\epsilon_j-\epsilon_r-\epsilon_s}$ with $\epsilon_i$ being the energy of orbital $i$ Nelec for the number of electrons, which equals the number of occupied spin orbitalsNorb for the total number of spin orbitalsNvirt = Norb - Nelec for the number of virtual (unoccupied) orbitalsorbitals are ordered by energy, thus occupied first, then virtualsAnd orbital indices follow the convention:a,b,c,d running over occupied orbitals onlyr,s,t,u running over virtual orbitals onlyi,j,k,l running over all orbitals
Optimizing for loops by using Einstein summations instead
python;numpy