fix: fill ERI in physicist notation in C

Author: Aayush Gupta

gbasis/integrals/libcint.py

@@ -1326,12 +1326,12 @@ def electron_repulsion(self):
             Electron repulsion integral array.
 
         """
-        out = np.zeros((self.nbfn, self.nbfn, self.nbfn, self.nbfn), dtype=c_double, order='F')
+        out = np.zeros((self.nbfn, self.nbfn, self.nbfn, self.nbfn), dtype=c_double)
         libcint_bindings.eri_shellloop(
             out, self.natm, self.atm, self.nbas,
             self.bas, self.env, self._offs, self.nbfn
         )
-        return out.transpose(0, 2, 1, 3)
+        return out
 
     def electron_repulsion_integral(self, notation="physicist", transform=None):
         r"""

gbasis/integrals/src/libcint_wrap.c

@@ -261,15 +261,14 @@ eri_shellloop(PyObject *self, PyObject *args)
                                 for (int s = 0; s < s_off; s++) {
                                     double val = buf[p + p_off*(q + q_off*(r + r_off*s))];
                                     int i = ipos+p, j = jpos+q, k = kpos+r, l = lpos+s;
-                                    out[i*nbfn*nbfn*nbfn + j*nbfn*nbfn + k*nbfn + l] = val;
-                                    out[i*nbfn*nbfn*nbfn + j*nbfn*nbfn + l*nbfn + k] = val;
-                                    out[j*nbfn*nbfn*nbfn + i*nbfn*nbfn + k*nbfn + l] = val;
-                                    out[j*nbfn*nbfn*nbfn + i*nbfn*nbfn + l*nbfn + k] = val;
-                                    out[k*nbfn*nbfn*nbfn + l*nbfn*nbfn + i*nbfn + j] = val;
-                                    out[k*nbfn*nbfn*nbfn + l*nbfn*nbfn + j*nbfn + i] = val;
-                                    out[l*nbfn*nbfn*nbfn + k*nbfn*nbfn + i*nbfn + j] = val;
-                                    out[l*nbfn*nbfn*nbfn + k*nbfn*nbfn + j*nbfn + i] = val;
-                                }
+                                    out[i*nbfn*nbfn*nbfn + k*nbfn*nbfn + j*nbfn + l] = val;
+                                    out[i*nbfn*nbfn*nbfn + l*nbfn*nbfn + j*nbfn + k] = val;
+                                    out[j*nbfn*nbfn*nbfn + k*nbfn*nbfn + i*nbfn + l] = val;
+                                    out[j*nbfn*nbfn*nbfn + l*nbfn*nbfn + i*nbfn + k] = val;
+                                    out[k*nbfn*nbfn*nbfn + i*nbfn*nbfn + l*nbfn + j] = val;
+                                    out[k*nbfn*nbfn*nbfn + j*nbfn*nbfn + l*nbfn + i] = val;
+                                    out[l*nbfn*nbfn*nbfn + i*nbfn*nbfn + k*nbfn + j] = val;
+                                    out[l*nbfn*nbfn*nbfn + j*nbfn*nbfn + k*nbfn + i] = val;                                }
                             }
                         }
                     }