Standardize Libcint implementation to include IOData conventions

.github/workflows/pytest.yaml

@@ -36,6 +36,10 @@ jobs:
           # Need editable mode in order to include the test files
           pip install -e .
 
+      - name: Install Libcint library
+        run: |
+          USE_LIBCINT=1 ./tools/install_libcint.sh
+
       - name: Run pytest
         uses: pavelzw/pytest-action@v2
         with:

.pre-commit-config.yaml

@@ -1,5 +1,6 @@
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
+  exclude: "notebooks/|tests\.fchk$"
   rev: v4.5.0
   hooks:
     - id: check-added-large-files                   # prevents giant files from being commit

gbasis/contractions.py

@@ -495,10 +495,13 @@ def angmom_components_sph(self):
         ['s2', 's1', 'c0', 'c1', 'c2'].
 
         """
-        return tuple(
-            ["s{}".format(m) for m in range(self.angmom, 0, -1)]
-            + ["c{}".format(m) for m in range(self.angmom + 1)]
-        )
+        if self.angmom == 1:
+            return tuple(["c1", "s1", "c0"])
+        else:
+            return tuple(
+                ["s{}".format(m) for m in range(self.angmom, 0, -1)]
+                + ["c{}".format(m) for m in range(self.angmom + 1)]
+            )
 
     @property
     def norm_prim_cart(self):

gbasis/integrals/libcint.py

@@ -458,6 +458,17 @@ def __init__(self, basis, atnums, atcoords, coord_type="spherical"):
         offs = np.asarray(offs, dtype=c_int)
         atm_offs = np.cumsum(atm_offs)
 
+        # Get permutation vector for ordering convention
+        permutations = []
+        for shell in basis:
+            if hasattr(shell, "permutation_libcint"):
+                permutation = shell.permutation_libcint()
+            else:
+                permutation = list(range(num_angmom(shell)))
+            for _ in range(shell.num_seg_cont):
+                perm_off = len(permutations)
+                permutations.extend(p + perm_off for p in permutation)
+
         # Allocate and fill C input arrays
         ienv = 20
         atm = np.zeros((natm, 6), dtype=c_int)
@@ -536,9 +547,10 @@ def __init__(self, basis, atnums, atcoords, coord_type="spherical"):
         self.atcoords = atcoords.copy()
         self._atm_offs = atm_offs
 
-        # Save basis function offsets
+        # Save basis function offsets and ordering permutation
         self._offs = offs
         self._max_off = max(offs)
+        self._permutations = permutations
 
         # Set inverse sqrt of overlap integral (temporarily, for __init__)
         self._ovlp_minhalf = np.ones(nbfn)
@@ -740,6 +752,9 @@ def int1e(notation="physicist", origin=None, inv_origin=None):
             if constant is not None:
                 out *= constant
 
+            # Apply permutation
+            out = out[self._permutations, :][:, self._permutations]
+
             # Return normalized integrals
             if self.coord_type == "cartesian":
                 return np.einsum(norm_einsum, self._ovlp_minhalf, self._ovlp_minhalf, out)
@@ -941,6 +956,12 @@ def int2e(notation="physicist", origin=None, inv_origin=None):
             if physicist:
                 out = out.transpose(0, 2, 1, 3)
 
+            # Apply permutation
+            out = out[self._permutations]
+            out = out[:, self._permutations]
+            out = out[:, :, self._permutations]
+            out = out[:, :, :, self._permutations]
+
             # Return normalized integrals
             if self.coord_type == "cartesian":
                 return np.einsum(

gbasis/wrappers.py

@@ -3,7 +3,28 @@
 import numpy as np
 
 
+CONVENTIONS_LIBCINT = {
+    (5, 'p'): ['s5','s4', 's3', 's2', 's1', 'c0', 'c1', 'c2', 'c3', 'c4','c5'],
+    (4, 'p'): ['s4','s3', 's2', 's1', 'c0', 'c1', 'c2', 'c3', 'c4'],
+    (3, 'p'): ['s3', 's2', 's1', 'c0', 'c1', 'c2', 'c3'],
+    (2, 'p'): ['s2', 's1', 'c0', 'c1', 'c2'],
+    (0, 'c'): ['1'],
+    (1, 'c'): ['x', 'y', 'z'],
+    (2, 'c'): ['xx', 'xy', 'xz', 'yy', 'yz', 'zz'],
+    (3, 'c'): ['xxx', 'xxy', 'xxz', 'xyy', 'xyz', 'xzz', 'yyy', 'yyz', 'yzz', 'zzz'],
+    (4, 'c'): ['xxxx', 'xxxy', 'xxxz', 'xxyy', 'xxyz', 'xxzz', 'xyyy', 'xyyz', 'xyzz',
+               'xzzz', 'yyyy', 'yyyz', 'yyzz', 'yzzz', 'zzzz'],
+    (5, 'c'): ['xxxxx', 'xxxxy', 'xxxxz', 'xxxyy', 'xxxyz', 'xxxzz', 'xxyyy', 'xxyyz', 'xxyzz',
+               'xxzzz', 'xyyyy', 'xyyyz', 'xyyzz', 'xyzzz', 'xzzzz', 'yyyyy', 'yyyyz', 'yyyzz',
+               'yyzzz', 'yzzzz', 'zzzzz'],
+    (6, 'c'): ['xxxxxx', 'xxxxxy', 'xxxxxz', 'xxxxyy', 'xxxxyz', 'xxxxzz', 'xxxyyy',
+               'xxxyyz', 'xxxyzz', 'xxxzzz', 'xxyyyy', 'xxyyyz', 'xxyyzz', 'xxyzzz',
+               'xxzzzz', 'xyyyyy', 'xyyyyz', 'xyyyzz', 'xyyzzz', 'xyzzzz', 'xzzzzz',
+               'yyyyyy', 'yyyyyz', 'yyyyzz', 'yyyzzz', 'yyzzzz', 'yzzzzz', 'zzzzzz'],
+}
+
 def from_iodata(mol, tol=1e-20, overlap=False):
+
     """Return basis set stored within the `IOData` instance in `iodata`.
 
     Parameters
@@ -113,6 +134,40 @@ def angmom_components_sph(self):
 
             return tuple(sph_conventions[self.angmom])
 
+        def permutation_libcint(self):
+            """Returns the permutation of the stored convention to the one used in libcint for
+            the given angular momentum
+
+            Conventions supported
+            cartesian: string of X, Y, and Z characters,such that :math:`a_x` `X` characters,
+            :math:`a_y` `Y` characters, and :math:`a_z` `Z`, characters appear. F
+            For example, `(2, 1, 1)` corresponds to `XXYZ`
+
+            Spherical: Strings of the form `c{m}` and `s{m}`.
+
+            Rerturns:
+            --------
+            permutation: An integer array that permutes basis function from IOData stored convention
+                         to LibCint.
+
+            """
+
+            angmom = self.angmom
+            if self.coord_type == "cartesian":
+                coord_type = 'c'
+                iodata_shell_convention = cart_conventions[angmom]
+            elif self.coord_type == 'spherical':
+                coord_type = 'p'
+                iodata_shell_convention = sph_conventions[angmom]
+
+            # Get libcint convention
+            libcint_convention = CONVENTIONS_LIBCINT[(angmom, coord_type)]
+            # Generate permutation
+            permutation = [libcint_convention.index(conv1) for conv1 in iodata_shell_convention]
+
+            return permutation
+
+
     if molbasis.primitive_normalization != "L2":  # pragma: no cover
         raise ValueError(
             "Only L2 normalization scheme is supported in `gbasis`. Given `IOData` instance uses "