From 02919109627e19e4af47670fc0275e3e288afb6a Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 22 Nov 2023 12:23:47 -0500
Subject: [PATCH 01/17] Update make_contractions

---
 gbasis/parsers.py | 28 ++++++++++++++++++++++++++--
 1 file changed, 26 insertions(+), 2 deletions(-)

diff --git a/gbasis/parsers.py b/gbasis/parsers.py
index 11d0e16f..17727d28 100644
--- a/gbasis/parsers.py
+++ b/gbasis/parsers.py
@@ -165,7 +165,7 @@ def parse_gbs(gbs_basis_file):
     return output
 
 
-def make_contractions(basis_dict, atoms, coords):
+def make_contractions(basis_dict, atoms, coords, coord_types='spherical'):
     """Return the contractions that correspond to the given atoms for the given basis.
 
     Parameters
@@ -176,6 +176,13 @@ def make_contractions(basis_dict, atoms, coords):
         Atoms at which the contractions are centered.
     coords : np.ndarray(N, 3)
         Coordinates of each atom.
+    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
+        Types of the coordinate system for the contractions.
+        If "cartesian", then all of the contractions are treated as Cartesian contractions.
+        If "spherical", then all of the contractions are treated as spherical contractions.
+        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
+        coordinate type of each `GeneralizedContractionShell` instance.
+        Default value is "spherical".
 
     Returns
     -------
@@ -202,7 +209,24 @@ def make_contractions(basis_dict, atoms, coords):
         raise ValueError("Number of atoms must be equal to the number of rows in the coordinates.")
 
     basis = []
+    len_coord_types = len(coord_types)
     for atom, coord in zip(atoms, coords):
         for angmom, exps, coeffs in basis_dict[atom]:
-            basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps))
+            if type(coord_types) == str:
+                # if coord_types given as a single string, assign the specified type to all contractions for all atoms
+                if coord_types == "spherical":
+                    basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, 'p'))
+                elif coord_types == "cartesian":
+                    basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, 'c'))
+                else:
+                    raise ValueError("If coord_types is a string, it must be either 'spherical' or 'cartesian'.")
+            elif type(coord_types) == list:
+                # if coord_types given as a list, assign the specified type to each atom's contractions individually
+                if len_coord_types == sum([len(basis_dict[i]) for i in atoms]):
+                    basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, coord_types.pop(0)))
+                else:
+                    raise ValueError("If coord_types is a list, it must be the same length as the total number of contractions.")
+            else:
+                raise TypeError("coord_types must be a string or list of strings.")
+
     return tuple(basis)

From 49e1b18dd9e664d5cf9456d636318812e0366ce4 Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 22 Nov 2023 12:27:51 -0500
Subject: [PATCH 02/17] Update docs and coord_type checks in base*

---
 gbasis/base.py           | 12 ++++++------
 gbasis/base_four_symm.py | 22 ++++++++++------------
 gbasis/base_one.py       | 22 ++++++++++------------
 gbasis/base_two_asymm.py | 30 ++++++++++++++----------------
 gbasis/base_two_symm.py  | 22 ++++++++++------------
 5 files changed, 50 insertions(+), 58 deletions(-)

diff --git a/gbasis/base.py b/gbasis/base.py
index a5e01019..6a3ab65c 100644
--- a/gbasis/base.py
+++ b/gbasis/base.py
@@ -157,12 +157,12 @@ def construct_array_lincomb(self, *transform, coord_type, **kwargs):
             Transformation matrix that will be used for linearly combining the spherical
             contractions.
             Note that multiple instances may be needed to construct the array.
-        coord_type : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
-            Types of the coordinate system for the contractions.
-            If "cartesian", then all of the contractions are treated as Cartesian contractions.
-            If "spherical", then all of the contractions are treated as spherical contractions.
-            If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-            coordinate type of each `GeneralizedContractionShell` instance.
+        coord_type : list/tuple of str
+            Types of the coordinate system for each GeneralizedContractionShell.
+            Each entry must be one of "cartesian" or "spherical". If multiple
+            instances of GeneralizedContractionShell are given but only one string
+            ("cartesian" or "spherical") is provided in the list/tuple, all of the
+            contractions will be treated according to that string.
         kwargs : dict
             Other keyword arguments that will be used to construct the array.
 
diff --git a/gbasis/base_four_symm.py b/gbasis/base_four_symm.py
index 1310dce9..03f8fa1c 100644
--- a/gbasis/base_four_symm.py
+++ b/gbasis/base_four_symm.py
@@ -568,12 +568,12 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
             Transformation is applied to the left.
             Rows correspond to the linear combinations (i.e. MO) and the columns correspond to the
             contractions (i.e. AO).
-        coord_type : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
-            Types of the coordinate system for the contractions.
-            If "cartesian", then all of the contractions are treated as Cartesian contractions.
-            If "spherical", then all of the contractions are treated as spherical contractions.
-            If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-            coordinate type of each `GeneralizedContractionShell` instance.
+        coord_type : list/tuple of str
+            Types of the coordinate system for each GeneralizedContractionShell.
+            Each entry must be one of "cartesian" or "spherical". If multiple
+            instances of GeneralizedContractionShell are given but only one string
+            ("cartesian" or "spherical") is provided in the list/tuple, all of the
+            contractions will be treated according to that string.
         kwargs : dict
             Other keyword arguments that will be used to construct the array.
             These keyword arguments are passed directly to `construct_array_spherical`, which will
@@ -592,20 +592,18 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
         Raises
         ------
         TypeError
-            If `coord_type` is not one of "cartesian", "spherical", or a list/tuple of these
-            strings.
+            If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if coord_type == "cartesian":
+        if all(type == "cartesian" for type in coord_type):
             array = self.construct_array_cartesian(**kwargs)
-        elif coord_type == "spherical":
+        elif all(type == "spherical" for type in coord_type):
             array = self.construct_array_spherical(**kwargs)
         elif isinstance(coord_type, (list, tuple)):
             array = self.construct_array_mix(coord_type, **kwargs)
         else:
             raise TypeError(
-                "`coord_type` must be one of 'cartesian', 'spherical', or a list/tuple of these "
-                "strings."
+                "`coord_type` must be a list/tuple of the strings 'cartesian' or 'spherical'"
             )
         array = np.tensordot(transform, array, (1, 0))
         array = np.tensordot(transform, array, (1, 1))
diff --git a/gbasis/base_one.py b/gbasis/base_one.py
index e05d13fe..e8270627 100644
--- a/gbasis/base_one.py
+++ b/gbasis/base_one.py
@@ -257,12 +257,12 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
             Transformation is applied to the left.
             Rows correspond to the linear combinationes (i.e. MO) and the columns correspond to the
             contractions (i.e. AO).
-        coord_type : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
-            Types of the coordinate system for the contractions.
-            If "cartesian", then all of the contractions are treated as Cartesian contractions.
-            If "spherical", then all of the contractions are treated as spherical contractions.
-            If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-            coordinate type of each GeneralizedContractionShell instance.
+        coord_type : list/tuple of str
+            Types of the coordinate system for each GeneralizedContractionShell.
+            Each entry must be one of "cartesian" or "spherical". If multiple
+            instances of GeneralizedContractionShell are given but only one string
+            ("cartesian" or "spherical") is provided in the list/tuple, all of the
+            contractions will be treated according to that string.
         kwargs : dict
             Other keyword arguments that will be used to construct the array.
             These keyword arguments are passed directly to `construct_array_spherical`, which will
@@ -278,19 +278,17 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
         Raises
         ------
         TypeError
-            If `coord_type` is not one of "cartesian", "spherical", or a list/tuple of these
-            strings.
+            If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if coord_type == "cartesian":
+        if all(type == "cartesian" for type in coord_type):
             array = self.construct_array_cartesian(**kwargs)
-        elif coord_type == "spherical":
+        elif all(type == "spherical" for type in coord_type):
             array = self.construct_array_spherical(**kwargs)
         elif isinstance(coord_type, (list, tuple)):
             array = self.construct_array_mix(coord_type, **kwargs)
         else:
             raise TypeError(
-                "`coord_type` must be one of 'cartesian', 'spherical', or a list/tuple of these "
-                "strings."
+                "`coord_type` must be a list/tuple of the strings 'cartesian' or 'spherical'"
             )
         return np.tensordot(transform, array, (1, 0))
diff --git a/gbasis/base_two_asymm.py b/gbasis/base_two_asymm.py
index 1926bf20..f1a86261 100644
--- a/gbasis/base_two_asymm.py
+++ b/gbasis/base_two_asymm.py
@@ -384,18 +384,18 @@ def construct_array_lincomb(
             Array associated with the linear combinations of spherical Gaussians (LCAO's) associated
             with the second index.
             If None, then transformation is skipped.
-        coord_type_one : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
+        coord_type_one : list/tuple of string
             Types of the coordinate system for the contractions associated with the first index.
-            If "cartesian", then all of the contractions are treated as Cartesian contractions.
-            If "spherical", then all of the contractions are treated as spherical contractions.
-            If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-            coordinate type of each `GeneralizedContractionShell` instance.
-        coord_type_two : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
+            Each entry must be one of "cartesian" or "spherical". If multiple
+            instances of GeneralizedContractionShell are given but only one string
+            ("cartesian" or "spherical") is provided in the list/tuple, all of the
+            contractions will be treated according to that string.
+        coord_type_two : list/tuple of string
             Types of the coordinate system for the contractions associated with the second index.
-            If "cartesian", then all of the contractions are treated as Cartesian contractions.
-            If "spherical", then all of the contractions are treated as spherical contractions.
-            If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-            coordinate type of each `GeneralizedContractionShell` instance.
+            Each entry must be one of "cartesian" or "spherical". If multiple
+            instances of GeneralizedContractionShell are given but only one string
+            ("cartesian" or "spherical") is provided in the list/tuple, all of the
+            contractions will be treated according to that string.
         kwargs : dict
             Other keyword arguments that will be used to construct the array.
             These keyword arguments are passed directly to `construct_array_spherical`, which will
@@ -418,13 +418,12 @@ def construct_array_lincomb(
         Raises
         ------
         TypeError
-            If `coord_type_one` and `coord_type_two` are not one of "cartesian", "spherical", or a
-            list/tuple of these strings.
+            If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if coord_type_one == "cartesian" and coord_type_two == "cartesian":
+        if all(type_one == "cartesian" for type_one in coord_type_one) and all(type_two == "cartesian" for type_two in coord_type_two):
             array = self.construct_array_cartesian(**kwargs)
-        elif coord_type_one == "spherical" and coord_type_two == "spherical":
+        elif all(type_one == "spherical" for type_one in coord_type_one) and all(type_two == "spherical" for type_two in coord_type_two):
             array = self.construct_array_spherical(**kwargs)
         else:
             if coord_type_one in ["cartesian", "spherical"]:
@@ -436,8 +435,7 @@ def construct_array_lincomb(
                 and isinstance(coord_type_two, (list, tuple))
             ):
                 raise TypeError(
-                    "`coord_type` must be one of 'cartesian', 'spherical', or a list/tuple of these"
-                    " strings."
+                    "`coord_type` must be a list/tuple of the strings 'cartesian' or 'spherical'"
                 )
             array = self.construct_array_mix(coord_type_one, coord_type_two, **kwargs)
         if transform_one is not None:
diff --git a/gbasis/base_two_symm.py b/gbasis/base_two_symm.py
index 7809d632..617914cb 100644
--- a/gbasis/base_two_symm.py
+++ b/gbasis/base_two_symm.py
@@ -364,12 +364,12 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
             Transformation is applied to the left.
             Rows correspond to the linear combinations (i.e. MO) and the columns correspond to the
             contractions (i.e. AO).
-        coord_type : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
-            Types of the coordinate system for the contractions.
-            If "cartesian", then all of the contractions are treated as Cartesian contractions.
-            If "spherical", then all of the contractions are treated as spherical contractions.
-            If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-            coordinate type of each `GeneralizedContractionShell` instance.
+        coord_type : list/tuple of str
+            Types of the coordinate system for each GeneralizedContractionShell.
+            Each entry must be one of "cartesian" or "spherical". If multiple
+            instances of GeneralizedContractionShell are given but only one string
+            ("cartesian" or "spherical") is provided in the list/tuple, all of the
+            contractions will be treated according to that string.
         kwargs : dict
             Other keyword arguments that will be used to construct the array.
             These keyword arguments are passed directly to `construct_array_spherical`, which will
@@ -388,20 +388,18 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
         Raises
         ------
         TypeError
-            If `coord_type` is not one of "cartesian", "spherical", or a list/tuple of these
-            strings.
+            If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if coord_type == "cartesian":
+        if all(type == "cartesian" for type in coord_type):
             array = self.construct_array_cartesian(**kwargs)
-        elif coord_type == "spherical":
+        elif all(type == "spherical" for type in coord_type):
             array = self.construct_array_spherical(**kwargs)
         elif isinstance(coord_type, (list, tuple)):
             array = self.construct_array_mix(coord_type, **kwargs)
         else:
             raise TypeError(
-                "`coord_type` must be one of 'cartesian', 'spherical', or a list/tuple of these "
-                "strings."
+                "`coord_type` must be a list/tuple of the strings 'cartesian' or 'spherical'"
             )
         array = np.tensordot(transform, array, (1, 0))
         array = np.tensordot(transform, array, (1, 1))

From 92ba2dbab8444882d042ae4a353f85e333a9683f Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 22 Nov 2023 12:52:24 -0500
Subject: [PATCH 03/17] Remove coord_type args from evals/integrals

---
 gbasis/evals/density.py                       | 90 +++----------------
 gbasis/evals/electrostatic_potential.py       | 17 ++--
 gbasis/evals/eval.py                          | 15 ++--
 gbasis/evals/eval_deriv.py                    | 15 ++--
 gbasis/evals/stress_tensor.py                 | 40 +--------
 gbasis/integrals/angular_momentum.py          | 19 ++--
 gbasis/integrals/electron_repulsion.py        | 15 ++--
 gbasis/integrals/kinetic_energy.py            | 16 ++--
 gbasis/integrals/moment.py                    | 16 ++--
 gbasis/integrals/momentum.py                  | 16 ++--
 .../integrals/nuclear_electron_attraction.py  | 11 +--
 gbasis/integrals/overlap.py                   | 17 ++--
 gbasis/integrals/overlap_asymm.py             | 19 +---
 gbasis/integrals/point_charge.py              | 15 ++--
 14 files changed, 73 insertions(+), 248 deletions(-)

diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py
index f36e80c2..9afb34b1 100644
--- a/gbasis/evals/density.py
+++ b/gbasis/evals/density.py
@@ -62,7 +62,7 @@ def evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval):
     return np.sum(density, axis=0)
 
 
-def evaluate_density(one_density_matrix, basis, points, transform=None, coord_type="spherical"):
+def evaluate_density(one_density_matrix, basis, points, transform=None):
     r"""Return the density of the given basis set at the given points.
 
     Parameters
@@ -84,13 +84,6 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, coord_ty
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -98,7 +91,7 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, coord_ty
         Density evaluated at `N` grid points.
 
     """
-    orb_eval = evaluate_basis(basis, points, transform=transform, coord_type=coord_type)
+    orb_eval = evaluate_basis(basis, points, transform=transform)
     return evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval)
 
 
@@ -109,7 +102,6 @@ def evaluate_deriv_reduced_density_matrix(
     basis,
     points,
     transform=None,
-    coord_type="spherical",
 ):
     r"""Return the derivative of the first-order reduced density matrix at the given points.
 
@@ -156,13 +148,6 @@ def evaluate_deriv_reduced_density_matrix(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -171,13 +156,13 @@ def evaluate_deriv_reduced_density_matrix(
 
     """
     deriv_orb_eval_one = evaluate_deriv_basis(
-        basis, points, orders_one, transform=transform, coord_type=coord_type
+        basis, points, orders_one, transform=transform
     )
     if np.array_equal(orders_one, orders_two):
         deriv_orb_eval_two = deriv_orb_eval_one
     else:
         deriv_orb_eval_two = evaluate_deriv_basis(
-            basis, points, orders_two, transform=transform, coord_type=coord_type
+            basis, points, orders_two, transform=transform
         )
     density = one_density_matrix.dot(deriv_orb_eval_two)
     density *= deriv_orb_eval_one
@@ -186,7 +171,7 @@ def evaluate_deriv_reduced_density_matrix(
 
 
 def evaluate_deriv_density(
-    orders, one_density_matrix, basis, points, transform=None, coord_type="spherical"
+    orders, one_density_matrix, basis, points, transform=None
 ):
     r"""Return the derivative of density of the given transformed basis set at the given points.
 
@@ -211,13 +196,6 @@ def evaluate_deriv_density(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -250,14 +228,13 @@ def evaluate_deriv_density(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
                 output += factor * num_occurence * density
     return output
 
 
 def evaluate_density_gradient(
-    one_density_matrix, basis, points, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, transform=None
 ):
     r"""Return the gradient of the density evaluated at the given points.
 
@@ -280,13 +257,6 @@ def evaluate_density_gradient(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -302,7 +272,6 @@ def evaluate_density_gradient(
                 basis,
                 points,
                 transform=transform,
-                coord_type=coord_type,
             ),
             evaluate_deriv_density(
                 np.array([0, 1, 0]),
@@ -310,7 +279,6 @@ def evaluate_density_gradient(
                 basis,
                 points,
                 transform=transform,
-                coord_type=coord_type,
             ),
             evaluate_deriv_density(
                 np.array([0, 0, 1]),
@@ -318,14 +286,13 @@ def evaluate_density_gradient(
                 basis,
                 points,
                 transform=transform,
-                coord_type=coord_type,
             ),
         ]
     ).T
 
 
 def evaluate_density_laplacian(
-    one_density_matrix, basis, points, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, transform=None
 ):
     r"""Return the Laplacian of the density evaluated at the given points.
 
@@ -348,13 +315,6 @@ def evaluate_density_laplacian(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -368,7 +328,6 @@ def evaluate_density_laplacian(
         basis,
         points,
         transform=transform,
-        coord_type=coord_type,
     )
     output += evaluate_deriv_density(
         np.array([0, 2, 0]),
@@ -376,7 +335,6 @@ def evaluate_density_laplacian(
         basis,
         points,
         transform=transform,
-        coord_type=coord_type,
     )
     output += evaluate_deriv_density(
         np.array([0, 0, 2]),
@@ -384,13 +342,12 @@ def evaluate_density_laplacian(
         basis,
         points,
         transform=transform,
-        coord_type=coord_type,
     )
     return output
 
 
 def evaluate_density_hessian(
-    one_density_matrix, basis, points, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, transform=None
 ):
     r"""Return the Hessian of the density evaluated at the given points.
 
@@ -413,13 +370,6 @@ def evaluate_density_hessian(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -439,7 +389,6 @@ def evaluate_density_hessian(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
                 for orders_one in np.identity(3, dtype=int)
             ]
@@ -449,7 +398,7 @@ def evaluate_density_hessian(
 
 
 def evaluate_posdef_kinetic_energy_density(
-    one_density_matrix, basis, points, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, transform=None
 ):
     r"""Return evaluations of positive definite kinetic energy density at the given points.
 
@@ -484,13 +433,6 @@ def evaluate_posdef_kinetic_energy_density(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -508,14 +450,13 @@ def evaluate_posdef_kinetic_energy_density(
             basis,
             points,
             transform=transform,
-            coord_type=coord_type,
         )
     return 0.5 * output
 
 
 # TODO: test against a reference
 def evaluate_general_kinetic_energy_density(
-    one_density_matrix, basis, points, alpha, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, alpha, transform=None
 ):
     r"""Return evaluations of general form of the kinetic energy density at the given points.
 
@@ -546,13 +487,6 @@ def evaluate_general_kinetic_energy_density(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -570,10 +504,10 @@ def evaluate_general_kinetic_energy_density(
         raise TypeError("`alpha` must be an int or float.")
 
     general_kinetic_energy_density = evaluate_posdef_kinetic_energy_density(
-        one_density_matrix, basis, points, transform=transform, coord_type=coord_type
+        one_density_matrix, basis, points, transform=transform
     )
     if alpha != 0:
         general_kinetic_energy_density += alpha * evaluate_density_laplacian(
-            one_density_matrix, basis, points, transform=transform, coord_type=coord_type
+            one_density_matrix, basis, points, transform=transform
         )
     return general_kinetic_energy_density
diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py
index dea0f19e..10b72b0c 100644
--- a/gbasis/evals/electrostatic_potential.py
+++ b/gbasis/evals/electrostatic_potential.py
@@ -10,7 +10,6 @@ def electrostatic_potential(
     nuclear_coords,
     nuclear_charges,
     transform=None,
-    coord_type="spherical",
     threshold_dist=0.0,
 ):
     r"""Return the electrostatic potentials of the basis set in the Cartesian form.
@@ -40,12 +39,6 @@ def electrostatic_potential(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
     threshold_dist : {float, 0.0}
         Threshold for rejecting nuclei whose distances to the points are less than the provided
         value. i.e. nuclei that are closer to the point than the threshold are discarded when
@@ -100,13 +93,15 @@ def electrostatic_potential(
     if threshold_dist < 0:
         raise ValueError("`threshold_dist` must be greater than or equal to zero.")
 
-    if coord_type == "cartesian":
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
+    if all(type == "cartesian" for type in coord_type):
         if sum(cont.num_cart * cont.num_seg_cont for cont in basis) != one_density_matrix.shape[0]:
             raise ValueError(
                 "`one_density_matrix` does not have number of rows/columns that is equal to the "
                 "total number of Cartesian contractions (atomic orbitals)."
             )
-    elif coord_type == "spherical":
+    elif all(type == "spherical" for type in coord_type):
         if sum(cont.num_sph * cont.num_seg_cont for cont in basis) != one_density_matrix.shape[0]:
             raise ValueError(
                 "`one_density_matrix` does not have number of rows/columns that is equal to the "
@@ -128,10 +123,10 @@ def electrostatic_potential(
             )
     else:
         raise TypeError(
-            "`coord_type` must be 'spherical', 'cartesian', or a list/tuple of these strings."
+            "`coord_type` must be a list/tuple of the strings 'spherical' or 'cartesian'."
         )
     hartree_potential = point_charge_integral(
-        basis, points, -np.ones(points.shape[0]), transform=transform, coord_type=coord_type
+        basis, points, -np.ones(points.shape[0]), transform=transform
     )
     hartree_potential *= one_density_matrix[:, :, None]
     hartree_potential = np.sum(hartree_potential, axis=(0, 1))
diff --git a/gbasis/evals/eval.py b/gbasis/evals/eval.py
index 734b0abb..54107012 100644
--- a/gbasis/evals/eval.py
+++ b/gbasis/evals/eval.py
@@ -111,7 +111,7 @@ def construct_array_contraction(contractions, points):
         return output
 
 
-def evaluate_basis(basis, points, transform=None, coord_type="spherical"):
+def evaluate_basis(basis, points, transform=None):
     r"""Evaluate the basis set in the given coordinate system at the given points.
 
     Parameters
@@ -129,13 +129,6 @@ def evaluate_basis(basis, points, transform=None, coord_type="spherical"):
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -147,10 +140,12 @@ def evaluate_basis(basis, points, transform=None, coord_type="spherical"):
         `N` is the number of coordinates at which the contractions are evaluated.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
     if transform is not None:
         return Eval(basis).construct_array_lincomb(transform, coord_type, points=points)
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return Eval(basis).construct_array_cartesian(points=points)
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return Eval(basis).construct_array_spherical(points=points)
     return Eval(basis).construct_array_mix(coord_type, points=points)
diff --git a/gbasis/evals/eval_deriv.py b/gbasis/evals/eval_deriv.py
index caac9a4f..a9cb7e6d 100644
--- a/gbasis/evals/eval_deriv.py
+++ b/gbasis/evals/eval_deriv.py
@@ -126,7 +126,7 @@ def construct_array_contraction(contractions, points, orders):
         return output
 
 
-def evaluate_deriv_basis(basis, points, orders, transform=None, coord_type="spherical"):
+def evaluate_deriv_basis(basis, points, orders, transform=None):
     r"""Evaluate the derivative of the basis set in the given coordinate system at the given points.
 
     Parameters
@@ -149,13 +149,6 @@ def evaluate_deriv_basis(basis, points, orders, transform=None, coord_type="sphe
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each GeneralizedContractionShell instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -167,12 +160,14 @@ def evaluate_deriv_basis(basis, points, orders, transform=None, coord_type="sphe
         `N` is the number of coordinates at which the contractions are evaluated.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
     if transform is not None:
         return EvalDeriv(basis).construct_array_lincomb(
             transform, coord_type, points=points, orders=orders
         )
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return EvalDeriv(basis).construct_array_cartesian(points=points, orders=orders)
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return EvalDeriv(basis).construct_array_spherical(points=points, orders=orders)
     return EvalDeriv(basis).construct_array_mix(coord_type, points=points, orders=orders)
diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py
index 0973f95f..9bd14ff4 100644
--- a/gbasis/evals/stress_tensor.py
+++ b/gbasis/evals/stress_tensor.py
@@ -9,7 +9,7 @@
 
 # TODO: need to be tested against reference
 def evaluate_stress_tensor(
-    one_density_matrix, basis, points, alpha=1, beta=0, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, alpha=1, beta=0, transform=None
 ):
     r"""Return the stress tensor evaluated at the given coordinates.
 
@@ -71,13 +71,6 @@ def evaluate_stress_tensor(
     beta : {int, float}
         Second parameter of the stress tensor.
         Default value is 0.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -107,7 +100,6 @@ def evaluate_stress_tensor(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
             if alpha != 1:
                 output[i, j] += (1 - alpha) * evaluate_deriv_reduced_density_matrix(
@@ -117,7 +109,6 @@ def evaluate_stress_tensor(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
             if i == j and beta != 0:
                 output[i, j] -= (
@@ -128,7 +119,6 @@ def evaluate_stress_tensor(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                 )
             output[j, i] = output[i, j]
@@ -137,7 +127,7 @@ def evaluate_stress_tensor(
 
 # TODO: need to be tested against reference
 def evaluate_ehrenfest_force(
-    one_density_matrix, basis, points, alpha=1, beta=0, transform=None, coord_type="spherical"
+    one_density_matrix, basis, points, alpha=1, beta=0, transform=None
 ):
     r"""Return the Ehrenfest force.
 
@@ -197,13 +187,6 @@ def evaluate_ehrenfest_force(
     beta : {int, float}
         Second parameter of the stress tensor.
         Default value is 0.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -232,7 +215,6 @@ def evaluate_ehrenfest_force(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
             if alpha != 1:
                 output[i] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix(
@@ -242,7 +224,6 @@ def evaluate_ehrenfest_force(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
             if alpha != 0.5:
                 output[i] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix(
@@ -252,7 +233,6 @@ def evaluate_ehrenfest_force(
                     basis,
                     points,
                     transform=transform,
-                    coord_type=coord_type,
                 )
             if beta != 0:
                 output[i] += (
@@ -264,7 +244,6 @@ def evaluate_ehrenfest_force(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                 )
     return output.T
@@ -278,7 +257,6 @@ def evaluate_ehrenfest_hessian(
     alpha=1,
     beta=0,
     transform=None,
-    coord_type="spherical",
     symmetric=False,
 ):
     r"""Return the Ehrenfest Hessian.
@@ -349,13 +327,6 @@ def evaluate_ehrenfest_hessian(
     beta : {int, float}
         Second parameter of the stress tensor.
         Default value is 0.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
     symmetric : {True, False}
         Flag for symmetrizing the Hessian.
         If True, then the Hessian is symmetrized by averaging it with its transpose.
@@ -390,7 +361,6 @@ def evaluate_ehrenfest_hessian(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                     output[i, j] += alpha * evaluate_deriv_reduced_density_matrix(
                         2 * orders_one,
@@ -399,7 +369,6 @@ def evaluate_ehrenfest_hessian(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                 if alpha != 1:
                     output[i, j] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix(
@@ -409,7 +378,6 @@ def evaluate_ehrenfest_hessian(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                     output[i, j] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix(
                         2 * orders_one + orders_two,
@@ -418,7 +386,6 @@ def evaluate_ehrenfest_hessian(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                 if alpha != 0.5:
                     output[i, j] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix(
@@ -428,7 +395,6 @@ def evaluate_ehrenfest_hessian(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                     output[i, j] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix(
                         orders_one + orders_two,
@@ -437,7 +403,6 @@ def evaluate_ehrenfest_hessian(
                         basis,
                         points,
                         transform=transform,
-                        coord_type=coord_type,
                     )
                 if beta != 0:
                     output[i, j] += (
@@ -449,7 +414,6 @@ def evaluate_ehrenfest_hessian(
                             basis,
                             points,
                             transform=transform,
-                            coord_type=coord_type,
                         )
                     )
     if symmetric:
diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py
index c45c4fb1..daefdd8a 100644
--- a/gbasis/integrals/angular_momentum.py
+++ b/gbasis/integrals/angular_momentum.py
@@ -47,7 +47,7 @@ class AngularMomentumIntegral(BaseTwoIndexSymmetric):
         Return the integral over the angular momentum operators associated with the contraction in
         the given coordinate system.
         :math:`K_cont` is the total number of contractions within the given basis set.
-    construct_array_spherical_lincomb(self, transform) : np.ndarray(K_orbs, K_orbs, 3)
+    construct_array_lincomb(self, transform, coord_type) : np.ndarray(K_orbs, K_orbs, 3)
         Return the integral over the angular momentum operator associated with linear combinations
         of spherical Gaussians (linear combinations of atomic orbitals).
         :math:`K_orbs` is the number of basis functions produced after the linear combinations.
@@ -155,7 +155,7 @@ def construct_array_contraction(contractions_one, contractions_two):
         return -1j * np.transpose(output, (3, 2, 4, 1, 0))
 
 
-def angular_momentum_integral(basis, transform=None, coord_type="spherical"):
+def angular_momentum_integral(basis, transform=None):
     r"""Return the integral over :math:`hat{L}` of the given basis set.
 
     Parameters
@@ -168,13 +168,6 @@ def angular_momentum_integral(basis, transform=None, coord_type="spherical"):
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", "spherical", list/tuple of "cartesian" or "spherical}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default is "spherical".
 
     Returns
     -------
@@ -185,10 +178,12 @@ def angular_momentum_integral(basis, transform=None, coord_type="spherical"):
         Dimension 2 corresponds to the direction of the angular momentum :math:`(x, y, z)`.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
     if transform is not None:
         return AngularMomentumIntegral(basis).construct_array_lincomb(transform, coord_type)
-    if coord_type == "spherical":
-        return AngularMomentumIntegral(basis).construct_array_spherical()
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return AngularMomentumIntegral(basis).construct_array_cartesian()
+    if all(type == "spherical" for type in coord_type):
+        return AngularMomentumIntegral(basis).construct_array_spherical()
     return AngularMomentumIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/electron_repulsion.py b/gbasis/integrals/electron_repulsion.py
index 024e9189..d3098de3 100644
--- a/gbasis/integrals/electron_repulsion.py
+++ b/gbasis/integrals/electron_repulsion.py
@@ -202,7 +202,7 @@ def construct_array_contraction(cls, cont_one, cont_two, cont_three, cont_four):
 
 
 def electron_repulsion_integral(
-    basis, transform=None, coord_type="spherical", notation="physicist"
+    basis, transform=None, notation="physicist"
 ):
     """Return the electron repulsion integrals fo the given basis set.
 
@@ -216,13 +216,6 @@ def electron_repulsion_integral(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
     notation : {"physicist", "chemist"}
         Convention with which the integrals are ordered.
         Default is Physicists' notation.
@@ -245,11 +238,13 @@ def electron_repulsion_integral(
     if notation not in ["physicist", "chemist"]:
         raise ValueError("`notation` must be one of 'physicist' or 'chemist'")
 
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
     if transform is not None:
         array = ElectronRepulsionIntegral(basis).construct_array_lincomb(transform, coord_type)
-    elif coord_type == "cartesian":
+    elif all(type == "cartesian" for type in coord_type):
         array = ElectronRepulsionIntegral(basis).construct_array_cartesian()
-    elif coord_type == "spherical":
+    elif all(type == "spherical" for type in coord_type):
         array = ElectronRepulsionIntegral(basis).construct_array_spherical()
     else:
         array = ElectronRepulsionIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py
index fac44d5b..8ee4b8a3 100644
--- a/gbasis/integrals/kinetic_energy.py
+++ b/gbasis/integrals/kinetic_energy.py
@@ -42,7 +42,7 @@ class KineticEnergyIntegral(BaseTwoIndexSymmetric):
         Return the kinetic energy integrals associated with the contraction in the given coordinate
         system.
         `K_cont` is the total number of contractions within the given basis set.
-    construct_array_spherical_lincomb(self, transform) : np.ndarray(K_orbs, K_orbs)
+    construct_array_lincomb(self, transform, coord_type) : np.ndarray(K_orbs, K_orbs)
         Return the kinetic energy integral associated with linear combinations of spherical
         Gaussians (linear combinations of atomic orbitals).
         `K_orbs` is the number of basis functions produced after the linear combinations.
@@ -118,7 +118,7 @@ def construct_array_contraction(contractions_one, contractions_two):
         return -0.5 * np.sum(output, axis=0)
 
 
-def kinetic_energy_integral(basis, transform=None, coord_type="spherical"):
+def kinetic_energy_integral(basis, transform=None):
     """Return kinetic energy integral of the given basis set.
 
     Parameters
@@ -131,13 +131,6 @@ def kinetic_energy_integral(basis, transform=None, coord_type="spherical"):
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -149,11 +142,12 @@ def kinetic_energy_integral(basis, transform=None, coord_type="spherical"):
         `K_orbs` is the total number of basis functions in the basis set.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return KineticEnergyIntegral(basis).construct_array_lincomb(transform, coord_type)
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return KineticEnergyIntegral(basis).construct_array_cartesian()
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return KineticEnergyIntegral(basis).construct_array_spherical()
     return KineticEnergyIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py
index 9c0a8c19..6f1291fb 100644
--- a/gbasis/integrals/moment.py
+++ b/gbasis/integrals/moment.py
@@ -40,7 +40,7 @@ class Moment(BaseTwoIndexSymmetric):
         Return the moment integrals associated with all of the contraction in the given coordinate
         system.
         `K_cont` is the total number of contractions within the given basis set.
-    construct_array_lincomb(self, transform) : np.ndarray(K_orbs, K_orbs)
+    construct_array_lincomb(self, transform, coord_type, **kwargs) : np.ndarray(K_orbs, K_orbs)
         Return the moment integrals associated with the linear combinations of contractions in the
         given coordinate system.
         `K_orbs` is the number of basis functions produced after the linear combinations.
@@ -157,7 +157,7 @@ class if the first four indices correspond to the segmented contraction and the
         return np.transpose(output, (1, 2, 3, 4, 0))
 
 
-def moment_integral(basis, moment_coord, moment_orders, transform=None, coord_type="spherical"):
+def moment_integral(basis, moment_coord, moment_orders, transform=None):
     """Return moment integral of the given basis set.
 
     Parameters
@@ -176,13 +176,6 @@ def moment_integral(basis, moment_coord, moment_orders, transform=None, coord_ty
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -201,16 +194,17 @@ def moment_integral(basis, moment_coord, moment_orders, transform=None, coord_ty
     order.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return Moment(basis).construct_array_lincomb(
             transform, coord_type, moment_coord=moment_coord, moment_orders=moment_orders
         )
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return Moment(basis).construct_array_cartesian(
             moment_coord=moment_coord, moment_orders=moment_orders
         )
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return Moment(basis).construct_array_spherical(
             moment_coord=moment_coord, moment_orders=moment_orders
         )
diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py
index 718e7ca6..f9755c80 100644
--- a/gbasis/integrals/momentum.py
+++ b/gbasis/integrals/momentum.py
@@ -44,7 +44,7 @@ class MomentumIntegral(BaseTwoIndexSymmetric):
         Return the integral over the momentum operator associated with the contraction in the given
         coordinate system.
         `K_cont` is the total number of contractions within the given basis set.
-    construct_array_spherical_lincomb(self, transform) : np.ndarray(K_orbs, K_orbs, 3)
+    construct_array_lincomb(self, transform, coord_type) : np.ndarray(K_orbs, K_orbs, 3)
         Return the integral over the momentum operator associated with linear combinations of
         spherical Gaussians (linear combinations of atomic orbitals).
         `K_orbs` is the number of basis functions produced after the linear combinations.
@@ -111,7 +111,7 @@ def construct_array_contraction(contractions_one, contractions_two):
         return -1j * np.transpose(output, (1, 2, 3, 4, 0))
 
 
-def momentum_integral(basis, transform=None, coord_type="spherical"):
+def momentum_integral(basis, transform=None):
     """Return integral over momentum operator of the given basis set.
 
     Parameters
@@ -124,13 +124,6 @@ def momentum_integral(basis, transform=None, coord_type="spherical"):
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -140,11 +133,12 @@ def momentum_integral(basis, transform=None, coord_type="spherical"):
         number of basis functions in the basis set.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return MomentumIntegral(basis).construct_array_lincomb(transform, coord_type)
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return MomentumIntegral(basis).construct_array_cartesian()
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return MomentumIntegral(basis).construct_array_spherical()
     return MomentumIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py
index 62ee0b91..355c2328 100644
--- a/gbasis/integrals/nuclear_electron_attraction.py
+++ b/gbasis/integrals/nuclear_electron_attraction.py
@@ -4,7 +4,7 @@
 
 
 def nuclear_electron_attraction_integral(
-    basis, nuclear_coords, nuclear_charges, transform=None, coord_type="spherical"
+    basis, nuclear_coords, nuclear_charges, transform=None
 ):
     """Return the nuclear electron attraction integrals of the basis set in the Cartesian form.
 
@@ -22,13 +22,6 @@ def nuclear_electron_attraction_integral(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -40,7 +33,7 @@ def nuclear_electron_attraction_integral(
     """
     return np.sum(
         point_charge_integral(
-            basis, nuclear_coords, nuclear_charges, transform=transform, coord_type=coord_type
+            basis, nuclear_coords, nuclear_charges, transform=transform
         ),
         axis=2,
     )
diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py
index 2f8c6c82..10eae345 100644
--- a/gbasis/integrals/overlap.py
+++ b/gbasis/integrals/overlap.py
@@ -40,7 +40,7 @@ class Overlap(BaseTwoIndexSymmetric):
     construct_array_mix(self, coord_types, **kwargs) : np.ndarray(K_cont, K_cont)
         Return the overlap integrals associated with the contraction in the given coordinate system.
         `K_cont` is the total number of contractions within the given basis set.
-    construct_array_lincomb(self, transform) : np.ndarray(K_orbs, K_orbs)
+    construct_array_lincomb(self, transform, coord_type) : np.ndarray(K_orbs, K_orbs)
         Return the overlap integrals associated with the linear combinations of contractions in the
         given coordinate system.
         `K_orbs` is the number of basis functions produced after the linear combinations.
@@ -107,7 +107,7 @@ def construct_array_contraction(contractions_one, contractions_two):
         )[0]
 
 
-def overlap_integral(basis, transform=None, coord_type="spherical"):
+def overlap_integral(basis, transform=None):
     """Return overlap integral of the given basis set.
 
     Parameters
@@ -120,13 +120,6 @@ def overlap_integral(basis, transform=None, coord_type="spherical"):
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -136,10 +129,12 @@ def overlap_integral(basis, transform=None, coord_type="spherical"):
         number of basis functions in the basis set.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
     if transform is not None:
         return Overlap(basis).construct_array_lincomb(transform, coord_type)
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return Overlap(basis).construct_array_cartesian()
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return Overlap(basis).construct_array_spherical()
     return Overlap(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/overlap_asymm.py b/gbasis/integrals/overlap_asymm.py
index 5e310054..4f0858ee 100644
--- a/gbasis/integrals/overlap_asymm.py
+++ b/gbasis/integrals/overlap_asymm.py
@@ -66,8 +66,6 @@ def overlap_integral_asymmetric(
     basis_two,
     transform_one=None,
     transform_two=None,
-    coord_type_one="spherical",
-    coord_type_two="spherical",
 ):
     """Return overlap integrals between two basis sets.
 
@@ -89,20 +87,6 @@ def overlap_integral_asymmetric(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type_one : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions in `basis_one`.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
-    coord_type_two : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions in `basis_two`.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -114,6 +98,9 @@ def overlap_integral_asymmetric(
         `K_orbs_2` is the number of basis functions in the `basis_two`.
 
     """
+    coord_type_one = [type for type in [shell.coord_type for shell in basis_one]]
+    coord_type_two = [type for type in [shell.coord_type for shell in basis_two]]
+
     return OverlapAsymmetric(basis_one, basis_two).construct_array_lincomb(
         transform_one, transform_two, coord_type_one, coord_type_two
     )
diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py
index af44e2a8..ab3d1857 100644
--- a/gbasis/integrals/point_charge.py
+++ b/gbasis/integrals/point_charge.py
@@ -266,7 +266,7 @@ def construct_array_contraction(
 
 
 def point_charge_integral(
-    basis, points_coords, points_charge, transform=None, coord_type="spherical"
+    basis, points_coords, points_charge, transform=None
 ):
     r"""Return the point-charge interaction integrals of basis set in the given coordinate systems.
 
@@ -286,13 +286,6 @@ def point_charge_integral(
         Transformation is applied to the left, i.e. the sum is over the index 1 of `transform`
         and index 0 of the array for contractions.
         Default is no transformation.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `GeneralizedContractionShell` instance.
-        Default value is "spherical".
 
     Returns
     -------
@@ -304,15 +297,17 @@ def point_charge_integral(
         `N` is the number of coordinates at which the contractions are evaluated.
 
     """
+    coord_type = [type for type in [shell.coord_type for shell in basis]]
+
     if transform is not None:
         return PointChargeIntegral(basis).construct_array_lincomb(
             transform, coord_type, points_coords=points_coords, points_charge=points_charge
         )
-    if coord_type == "cartesian":
+    if all(type == "cartesian" for type in coord_type):
         return PointChargeIntegral(basis).construct_array_cartesian(
             points_coords=points_coords, points_charge=points_charge
         )
-    if coord_type == "spherical":
+    if all(type == "spherical" for type in coord_type):
         return PointChargeIntegral(basis).construct_array_spherical(
             points_coords=points_coords, points_charge=points_charge
         )

From fa125cc2cb094ae4b8a0829b7e51a7e807fce269 Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 22 Nov 2023 13:05:15 -0500
Subject: [PATCH 04/17] Add coord_type attr to GeneralizedContractionShell

---
 gbasis/contractions.py | 53 +++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 52 insertions(+), 1 deletion(-)

diff --git a/gbasis/contractions.py b/gbasis/contractions.py
index eca3111d..06887495 100644
--- a/gbasis/contractions.py
+++ b/gbasis/contractions.py
@@ -91,6 +91,8 @@ class GeneralizedContractionShell:
         Contraction coefficients, :math:`\{d_{ij}\}`, of the primitives.
         First axis corresponds to the primitive and the second axis corresponds to the segmented
         contraction shell.
+    coord_type : str
+        Type of the coordinate system used to specify the contractions.
     norm_cont : np.ndarray(M, L)
         Normalization constants of the Cartesian contractions of different angular momentum
         components and segmented contraction shells.
@@ -114,7 +116,7 @@ class GeneralizedContractionShell:
 
     """
 
-    def __init__(self, angmom, coord, coeffs, exps):
+    def __init__(self, angmom, coord, coeffs, exps, coord_type='p'):
         r"""Initialize a GeneralizedContractionShell instance.
 
         Parameters
@@ -135,6 +137,9 @@ def __init__(self, angmom, coord, coeffs, exps):
             dimension.
         exps : np.ndarray(K,)
             Exponents of the primitives, :math:`\{\alpha_i\}_{i=1}^K`.
+        coord_type : str
+            Coordinate type of the contraction.
+            Default is 'p' (spherical).
 
         """
         self.angmom = angmom
@@ -142,6 +147,7 @@ def __init__(self, angmom, coord, coeffs, exps):
         self.coeffs = coeffs
         self.exps = exps
         self.assign_norm_cont()
+        self.coord_type = coord_type
 
     @property
     def coord(self):
@@ -479,3 +485,48 @@ def assign_norm_cont(self):
 
         self.norm_cont = np.einsum("ijij->ij", Overlap.construct_array_contraction(self, self))
         self.norm_cont **= -0.5
+
+    @property
+    def coord_type(self):
+        """Return the coordinate type.
+
+        Returns
+        -------
+        coord_type : str
+            Coordinate type of the contraction.
+
+        """
+        return self._coord_type
+
+    @coord_type.setter
+    def coord_type(self, coord_type):
+        """Set the coordinate type.
+
+        Parameters
+        ----------
+        coord_type : str
+            Coordinate type of the contraction.
+
+        Raises
+        ------
+        TypeError
+            If `coord_type` is not a string.
+        ValueError
+            If `coord_type` is not one of the following:
+                - 'c' (cartesian)
+                - 'cartesian'
+                - 'p' (spherical)
+                - 'spherical'
+        """
+
+
+        if not isinstance(coord_type, str):
+            raise TypeError("Coordinate type must be given as a string.")
+        if coord_type not in ["c", "cartesian", "p", "spherical"]:
+            raise ValueError("`coord_type` is incorrectly specified. It must be either 'c' "
+                             "or 'cartesian' for Cartesian coordinates, or 'p' or 'spherical' "
+                             "for spherical coordinates.")
+        self._coord_type = {"c": "cartesian",
+                            "cartesian": "cartesian",
+                            "spherical": "spherical",
+                            "p": "spherical"}[coord_type]

From 71cda1e1b12eb1c68e00f185dc5991441fd7826f Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 22 Nov 2023 13:06:46 -0500
Subject: [PATCH 05/17] Negate from_iodata return of coord_types

---
 gbasis/wrappers.py | 13 ++-----------
 1 file changed, 2 insertions(+), 11 deletions(-)

diff --git a/gbasis/wrappers.py b/gbasis/wrappers.py
index da274ad2..b75cdffc 100644
--- a/gbasis/wrappers.py
+++ b/gbasis/wrappers.py
@@ -16,11 +16,6 @@ def from_iodata(mol):
     basis : tuple of gbasis.contraciton.GeneralizedContractionShell
         Basis set object used within the `gbasis` module.
         `GeneralizedContractionShell` corresponds to the `Shell` object within `iodata.basis`.
-    coord_types : list of str
-        A list of strings, each on being either ``"cartesian"`` or
-        ``"spherical"``, which needs to be passed into integral and eval
-        functions. This part of the basis set is not stored in the ``basis``
-        return value. Users need to keep track of this in a separate variable.
 
     Raises
     ------
@@ -127,7 +122,6 @@ def angmom_components_sph(self):
         )
 
     basis = []
-    coord_types = []
     for shell in molbasis.shells:
         # Verify that this is not a generalized contraction.
         if shell.ncon != 1:
@@ -136,15 +130,12 @@ def angmom_components_sph(self):
         # NOTE: GeneralizedContractionShell only accepts angular momentum as an int.
         angmom = int(shell.angmoms[0])
 
-        # get type
-        coord_types.append({"c": "cartesian", "p": "spherical"}[shell.kinds[0]])
-
         # pylint: disable=E1136
         basis.append(
-            IODataShell(angmom, mol.atcoords[shell.icenter], shell.coeffs, shell.exponents)
+            IODataShell(angmom, mol.atcoords[shell.icenter], shell.coeffs, shell.exponents, shell.kinds[0])
         )
 
-    return basis, coord_types
+    return basis
 
 
 def from_pyscf(mol):

From a5670a98414109922036c63c21116e7828967f8b Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 22 Nov 2023 13:18:56 -0500
Subject: [PATCH 06/17] Adjust tests for coord_types API change

---
 tests/test_angular_momentum.py            | 18 +++---
 tests/test_base_four_symm.py              |  6 +-
 tests/test_base_one.py                    | 10 +--
 tests/test_base_two_symm.py               | 14 ++--
 tests/test_electron_repulsion.py          | 46 ++++++-------
 tests/test_electrostatic_potential.py     | 79 +++++++++--------------
 tests/test_eval.py                        | 57 +++++++++-------
 tests/test_eval_deriv.py                  | 45 +++++++------
 tests/test_kinetic_energy.py              | 30 ++++-----
 tests/test_moment.py                      | 11 ++--
 tests/test_momentum.py                    | 18 +++---
 tests/test_nuclear_electron_attraction.py | 32 ++++-----
 tests/test_overlap.py                     | 33 +++++-----
 tests/test_overlap_asymm.py               | 45 +++++++------
 tests/test_point_charge.py                | 16 ++---
 tests/test_wrappers.py                    | 15 +++--
 16 files changed, 235 insertions(+), 240 deletions(-)

diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py
index 1acf5412..9f73d81f 100644
--- a/tests/test_angular_momentum.py
+++ b/tests/test_angular_momentum.py
@@ -369,11 +369,11 @@ def test_angular_momentum_construct_array_contraction():
 def test_angular_momentum_integral_cartesian():
     """Test gbasis.integrals.angular_momentum.angular_momentum_integral_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     angular_momentum_integral_obj = AngularMomentumIntegral(basis)
     assert np.allclose(
         angular_momentum_integral_obj.construct_array_cartesian(),
-        angular_momentum_integral(basis, coord_type="cartesian"),
+        angular_momentum_integral(basis),
     )
 
 
@@ -381,11 +381,11 @@ def test_angular_momentum_integral_spherical():
     """Test gbasis.integrals.angular_momentum.angular_momentum_integral_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     angular_momentum_integral_obj = AngularMomentumIntegral(basis)
     assert np.allclose(
         angular_momentum_integral_obj.construct_array_spherical(),
-        angular_momentum_integral(basis, coord_type="spherical"),
+        angular_momentum_integral(basis),
     )
 
 
@@ -393,21 +393,21 @@ def test_angular_momentum_integral_mix():
     """Test gbasis.integrals.angular_momentum.angular_momentum_integral_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ['spherical'] * 8)
     angular_momentum_integral_obj = AngularMomentumIntegral(basis)
     assert np.allclose(
         angular_momentum_integral_obj.construct_array_mix(["spherical"] * 8),
-        angular_momentum_integral(basis, coord_type=["spherical"] * 8),
+        angular_momentum_integral(basis),
     )
 
 
 def test_angular_momentum_integral_lincomb():
     """Test gbasis.integrals.angular_momentum.angular_momentum_integral_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     angular_momentum_integral_obj = AngularMomentumIntegral(basis)
     transform = np.random.rand(14, 18)
     assert np.allclose(
-        angular_momentum_integral_obj.construct_array_lincomb(transform, "spherical"),
-        angular_momentum_integral(basis, transform, coord_type="spherical"),
+        angular_momentum_integral_obj.construct_array_lincomb(transform, ["spherical"]),
+        angular_momentum_integral(basis, transform),
     )
diff --git a/tests/test_base_four_symm.py b/tests/test_base_four_symm.py
index 9f49da11..381e18bd 100644
--- a/tests/test_base_four_symm.py
+++ b/tests/test_base_four_symm.py
@@ -676,7 +676,7 @@ def test_construct_array_lincomb():
     contractions.norm_cont = np.ones((1, 3))
     test = Test([contractions])
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "cartesian"),
+        test.construct_array_lincomb(orb_transform, ["cartesian"]),
         np.einsum(
             "ijkl,ai,bj,ck,dl->abcd",
             np.einsum("ijkl->lkji", np.arange(81).reshape(3, 3, 3, 3)) * 2,
@@ -687,7 +687,7 @@ def test_construct_array_lincomb():
         ),
     )
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical"),
+        test.construct_array_lincomb(orb_transform, ["spherical"]),
         np.einsum(
             "ijkl,ai,bj,ck,dl->abcd",
             np.einsum(
@@ -1002,7 +1002,7 @@ def test_construct_array_lincomb():
     )
 
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical"),
+        test.construct_array_lincomb(orb_transform, ["spherical"]),
         np.einsum(
             "ijkl,ai,bj,ck,dl->abcd",
             sph_array,
diff --git a/tests/test_base_one.py b/tests/test_base_one.py
index 0a644825..043166b7 100644
--- a/tests/test_base_one.py
+++ b/tests/test_base_one.py
@@ -277,26 +277,26 @@ def test_contruct_array_lincomb():
     )
     test = Test([contractions])
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "cartesian"),
+        test.construct_array_lincomb(orb_transform, ["cartesian"]),
         orb_transform.dot(np.arange(9).reshape(3, 3)) * 2,
     )
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical"),
+        test.construct_array_lincomb(orb_transform, ["spherical"]),
         orb_transform.dot(sph_transform).dot(np.arange(9).reshape(3, 3)) * 2,
     )
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical", a=3),
+        test.construct_array_lincomb(orb_transform, ["spherical"], a=3),
         orb_transform.dot(sph_transform).dot(np.arange(9).reshape(3, 3)) * 3,
     )
     with pytest.raises(TypeError):
         test.construct_array_lincomb(orb_transform, "bad")
     with pytest.raises(TypeError):
-        test.construct_array_lincomb(orb_transform, "spherical", bad_keyword=3)
+        test.construct_array_lincomb(orb_transform, ["spherical"], bad_keyword=3)
 
     orb_transform = np.random.rand(3, 6)
     test = Test([contractions, contractions])
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical"),
+        test.construct_array_lincomb(orb_transform, ["spherical"]),
         orb_transform.dot(
             np.vstack([sph_transform.dot(np.arange(9, dtype=float).reshape(3, 3)) * 2] * 2)
         ),
diff --git a/tests/test_base_two_symm.py b/tests/test_base_two_symm.py
index 0aee5381..a5add74d 100644
--- a/tests/test_base_two_symm.py
+++ b/tests/test_base_two_symm.py
@@ -507,11 +507,11 @@ def test_contruct_array_lincomb():
     contractions.norm_cont = np.ones((1, 3))
     test = Test([contractions])
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "cartesian"),
+        test.construct_array_lincomb(orb_transform, ["cartesian"]),
         (orb_transform.dot(np.arange(9).reshape(3, 3)).dot(orb_transform.T).T * 2),
     )
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical"),
+        test.construct_array_lincomb(orb_transform, ["spherical"]),
         (
             orb_transform.dot(sph_transform)
             .dot(np.arange(9).reshape(3, 3))
@@ -522,7 +522,7 @@ def test_contruct_array_lincomb():
         ),
     )
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical", a=3),
+        test.construct_array_lincomb(orb_transform, ["spherical"], a=3),
         (
             orb_transform.dot(sph_transform)
             .dot(np.arange(9).reshape(3, 3))
@@ -533,7 +533,7 @@ def test_contruct_array_lincomb():
         ),
     )
     with pytest.raises(TypeError):
-        test.construct_array_lincomb(orb_transform, "spherical", bad_keyword=3)
+        test.construct_array_lincomb(orb_transform, ["spherical"], bad_keyword=3)
     with pytest.raises(TypeError):
         test.construct_array_lincomb(orb_transform, "bad", keyword=3)
 
@@ -561,7 +561,7 @@ def test_contruct_array_lincomb():
     )
     orb_transform = np.random.rand(8, 8)
     assert np.allclose(
-        test.construct_array_lincomb(orb_transform, "spherical", a=4),
+        test.construct_array_lincomb(orb_transform, ["spherical"], a=4),
         np.swapaxes(
             np.tensordot(
                 orb_transform,
@@ -830,13 +830,13 @@ def construct_array_contraction(self, cont_one, cont_two, a=2):
         test_symm.construct_array_spherical(), test_asymm.construct_array_spherical()
     )
     assert np.allclose(
-        test_symm.construct_array_lincomb(sph_orb_transform, "spherical"),
+        test_symm.construct_array_lincomb(sph_orb_transform, ["spherical"]),
         test_asymm.construct_array_lincomb(
             sph_orb_transform, sph_orb_transform, "spherical", "spherical"
         ),
     )
     assert np.allclose(
-        test_symm.construct_array_lincomb(cart_orb_transform, "cartesian"),
+        test_symm.construct_array_lincomb(cart_orb_transform, ["cartesian"]),
         test_asymm.construct_array_lincomb(
             cart_orb_transform, cart_orb_transform, "cartesian", "cartesian"
         ),
diff --git a/tests/test_electron_repulsion.py b/tests/test_electron_repulsion.py
index 6077b817..c76c73b0 100644
--- a/tests/test_electron_repulsion.py
+++ b/tests/test_electron_repulsion.py
@@ -99,12 +99,12 @@ def test_electron_repulsion_cartesian_horton_sto6g_bec():
     """
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
     coords = np.array([[0, 0, 0], [1.0, 0, 0]])
-    basis = make_contractions(basis_dict, ["Be", "C"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["Be", "C"], coords, 'cartesian')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_elec_repulsion = np.load(find_datafile("data_horton_bec_cart_elec_repulsion.npy"))
     assert np.allclose(
-        horton_elec_repulsion, electron_repulsion_integral(basis, coord_type="cartesian")
+        horton_elec_repulsion, electron_repulsion_integral(basis)
     )
 
 
@@ -120,8 +120,8 @@ def test_electron_repulsion_cartesian_horton_custom_hhe():
     """
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     coords = np.array([[0, 0, 0], [0.8, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs[:, 0], i.exps) for i in basis[:8]]
+    basis = make_contractions(basis_dict, ["H", "He"], coords, 'cartesian')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs[:, 0], i.exps, i.coord_type) for i in basis[:8]]
     basis[0] = HortonContractions(
         basis[0].angmom, basis[0].coord, basis[0].coeffs[3:], basis[0].exps[3:]
     )
@@ -133,62 +133,62 @@ def test_electron_repulsion_cartesian_horton_custom_hhe():
 
     horton_elec_repulsion = np.load(find_datafile("data_horton_hhe_cart_elec_repulsion.npy"))
     assert np.allclose(
-        horton_elec_repulsion, electron_repulsion_integral(basis, coord_type="cartesian")
+        horton_elec_repulsion, electron_repulsion_integral(basis)
     )
 
 
 def test_electron_repulsion_cartesian():
     """Test gbasis.integrals.electron_repulsion.electron_repulsion_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]), 'cartesian')
 
     erep_obj = ElectronRepulsionIntegral(basis)
     assert np.allclose(
         erep_obj.construct_array_cartesian(),
-        electron_repulsion_integral(basis, notation="chemist", coord_type="cartesian"),
+        electron_repulsion_integral(basis, notation="chemist"),
     )
     assert np.allclose(
         np.einsum("ijkl->ikjl", erep_obj.construct_array_cartesian()),
-        electron_repulsion_integral(basis, notation="physicist", coord_type="cartesian"),
+        electron_repulsion_integral(basis, notation="physicist"),
     )
     with pytest.raises(ValueError):
-        electron_repulsion_integral(basis, notation="bad", coord_type="cartesian")
+        electron_repulsion_integral(basis, notation="bad")
 
 
 def test_electron_repulsion_spherical():
     """Test gbasis.integrals.electron_repulsion.electron_repulsion_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]), 'spherical')
 
     erep_obj = ElectronRepulsionIntegral(basis)
     assert np.allclose(
         erep_obj.construct_array_spherical(),
-        electron_repulsion_integral(basis, notation="chemist", coord_type="spherical"),
+        electron_repulsion_integral(basis, notation="chemist"),
     )
     assert np.allclose(
         np.einsum("ijkl->ikjl", erep_obj.construct_array_spherical()),
-        electron_repulsion_integral(basis, notation="physicist", coord_type="spherical"),
+        electron_repulsion_integral(basis, notation="physicist"),
     )
     with pytest.raises(ValueError):
-        electron_repulsion_integral(basis, notation="bad", coord_type="spherical")
+        electron_repulsion_integral(basis, notation="bad")
 
 
 def test_electron_repulsion_mix():
     """Test gbasis.integrals.electron_repulsion.electron_repulsion_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]), ["spherical"] * 3)
 
     erep_obj = ElectronRepulsionIntegral(basis)
     assert np.allclose(
         erep_obj.construct_array_mix(["spherical"] * 3),
-        electron_repulsion_integral(basis, notation="chemist", coord_type=["spherical"] * 3),
+        electron_repulsion_integral(basis, notation="chemist"),
     )
     assert np.allclose(
         np.einsum("ijkl->ikjl", erep_obj.construct_array_mix(["spherical"] * 3)),
-        electron_repulsion_integral(basis, notation="physicist", coord_type=["spherical"] * 3),
+        electron_repulsion_integral(basis, notation="physicist"),
     )
     with pytest.raises(ValueError):
-        electron_repulsion_integral(basis, notation="bad", coord_type=["spherical"] * 3)
+        electron_repulsion_integral(basis, notation="bad")
 
 
 def test_electron_repulsion_lincomb():
@@ -199,12 +199,12 @@ def test_electron_repulsion_lincomb():
     erep_obj = ElectronRepulsionIntegral(basis)
     transform = np.random.rand(3, 5)
     assert np.allclose(
-        erep_obj.construct_array_lincomb(transform, "spherical"),
-        electron_repulsion_integral(basis, transform, notation="chemist", coord_type="spherical"),
+        erep_obj.construct_array_lincomb(transform, ["spherical"]),
+        electron_repulsion_integral(basis, transform, notation="chemist"),
     )
     assert np.allclose(
-        np.einsum("ijkl->ikjl", erep_obj.construct_array_lincomb(transform, "spherical")),
-        electron_repulsion_integral(basis, transform, notation="physicist", coord_type="spherical"),
+        np.einsum("ijkl->ikjl", erep_obj.construct_array_lincomb(transform, ["spherical"])),
+        electron_repulsion_integral(basis, transform, notation="physicist"),
     )
     with pytest.raises(ValueError):
-        electron_repulsion_integral(basis, transform, notation="bad", coord_type="spherical")
+        electron_repulsion_integral(basis, transform, notation="bad")
diff --git a/tests/test_electrostatic_potential.py b/tests/test_electrostatic_potential.py
index a8a17532..2d77214e 100644
--- a/tests/test_electrostatic_potential.py
+++ b/tests/test_electrostatic_potential.py
@@ -14,144 +14,125 @@ def test_electrostatic_potential():
     """
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    cartesian_basis = make_contractions(basis_dict, ["H", "He"], coords, 'cartesian')
+    cartesian_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in cartesian_basis]
+    spherical_basis = make_contractions(basis_dict, ["H", "He"], coords, 'spherical')
+    spherical_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in spherical_basis]
+    mixed_basis = make_contractions(basis_dict, ["H", "He"], coords, ['spherical'] * 9)
+    mixed_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in mixed_basis]
 
     # check density_matrix type
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103).tolist(),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103).flatten(),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     # check nuclear_coords
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3).tolist(),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 4),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     # check nuclear charges
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]).reshape(1, 2),
-            coord_type="cartesian",
         )
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]).tolist(),
-            coord_type="cartesian",
         )
     # check density_matrix symmetry
     with pytest.raises(ValueError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.eye(103, 102),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     # check nuclear_coords and nuclear_charges shapes
     with pytest.raises(ValueError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(3, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     # check threshold_dist types
     with pytest.raises(TypeError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
             threshold_dist=None,
         )
     # check threshold_dist value
     with pytest.raises(ValueError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
             threshold_dist=-0.1,
         )
-    # check coord_types type
-    with pytest.raises(TypeError):
-        electrostatic_potential(
-            basis,
-            np.identity(103),
-            np.random.rand(10, 3),
-            np.random.rand(2, 3),
-            np.array([1, 2]),
-            coord_type="bad",
-        )
     with pytest.raises(ValueError):
         electrostatic_potential(
-            basis,
+            cartesian_basis,
             np.identity(88),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="cartesian",
         )
     with pytest.raises(ValueError):
         electrostatic_potential(
-            basis,
+            spherical_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type="spherical",
         )
     with pytest.raises(ValueError):
         electrostatic_potential(
-            basis,
+            mixed_basis,
             np.identity(103),
             np.random.rand(10, 3),
             np.random.rand(2, 3),
             np.array([1, 2]),
-            coord_type=["spherical"] * 9,
         )
 
 
@@ -165,8 +146,8 @@ def test_electrostatic_potential_cartesian():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], coords, 'cartesian')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     grid_1d = np.linspace(-2, 2, num=5)
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
@@ -175,7 +156,7 @@ def test_electrostatic_potential_cartesian():
     horton_nucattract = np.load(find_datafile("data_horton_hhe_cart_esp.npy"))
     assert np.allclose(
         electrostatic_potential(
-            basis, np.identity(103), grid_3d, coords, np.array([1, 2]), coord_type="cartesian"
+            basis, np.identity(103), grid_3d, coords, np.array([1, 2])
         ),
         horton_nucattract,
     )
@@ -191,8 +172,8 @@ def test_electrostatic_potential_spherical():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], coords, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     grid_1d = np.linspace(-2, 2, num=5)
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
@@ -201,7 +182,7 @@ def test_electrostatic_potential_spherical():
     horton_nucattract = np.load(find_datafile("data_horton_hhe_sph_esp.npy"))
     assert np.allclose(
         electrostatic_potential(
-            basis, np.identity(88), grid_3d, coords, np.array([1, 2]), coord_type="spherical"
+            basis, np.identity(88), grid_3d, coords, np.array([1, 2])
         ),
         horton_nucattract,
     )
@@ -217,8 +198,10 @@ def test_electrostatic_potential_mix():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    spherical_basis = make_contractions(basis_dict, ["H", "He"], coords, ['spherical'] * 9)
+    spherical_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in spherical_basis]
+    cartesian_basis = make_contractions(basis_dict, ["H", "He"], coords, ['cartesian'] * 9)
+    cartesian_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in cartesian_basis]
 
     grid_1d = np.linspace(-2, 2, num=5)
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
@@ -227,14 +210,14 @@ def test_electrostatic_potential_mix():
     horton_nucattract = np.load(find_datafile("data_horton_hhe_sph_esp.npy"))
     assert np.allclose(
         electrostatic_potential(
-            basis, np.identity(88), grid_3d, coords, np.array([1, 2]), coord_type=["spherical"] * 9
+            spherical_basis, np.identity(88), grid_3d, coords, np.array([1, 2])
         ),
         horton_nucattract,
     )
     horton_nucattract = np.load(find_datafile("data_horton_hhe_cart_esp.npy"))
     assert np.allclose(
         electrostatic_potential(
-            basis, np.identity(103), grid_3d, coords, np.array([1, 2]), coord_type=["cartesian"] * 9
+            cartesian_basis, np.identity(103), grid_3d, coords, np.array([1, 2])
         ),
         horton_nucattract,
     )
diff --git a/tests/test_eval.py b/tests/test_eval.py
index 638e2ec3..1765c5b4 100644
--- a/tests/test_eval.py
+++ b/tests/test_eval.py
@@ -56,11 +56,11 @@ def test_evaluate_construct_array_contraction():
 def test_evaluate_basis_cartesian():
     """Test gbasis.evals.eval.evaluate_basis_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), 'cartesian')
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(points=np.array([[0, 0, 0]])),
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type="cartesian"),
+        evaluate_basis(basis, np.array([[0, 0, 0]])),
     )
 
 
@@ -69,25 +69,25 @@ def test_evaluate_basis_spherical():
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
     # cartesian and spherical are the same for s orbital
-    basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), 'spherical')
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(points=np.array([[0, 0, 0]])),
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type="spherical"),
+        evaluate_basis(basis, np.array([[0, 0, 0]])),
     )
     # p orbitals are zero at center
-    basis = make_contractions(basis_dict, ["Li"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Li"], np.array([[0, 0, 0]]), 'spherical')
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(points=np.array([[0, 0, 0]])),
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type="spherical"),
+        evaluate_basis(basis, np.array([[0, 0, 0]])),
     )
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_spherical(points=np.array([[1, 1, 1]])),
-        evaluate_basis(basis, np.array([[1, 1, 1]]), coord_type="spherical"),
+        evaluate_basis(basis, np.array([[1, 1, 1]])),
     )
 
 
@@ -96,36 +96,44 @@ def test_evaluate_basis_mix():
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
     # cartesian and spherical are the same for s orbital
-    basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]))
+    spherical_basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), "spherical")
+    spherical_basis_list = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), ["spherical"])
     assert np.allclose(
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type="spherical"),
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type=["spherical"]),
+        evaluate_basis(spherical_basis, np.array([[0, 0, 0]])),
+        evaluate_basis(spherical_basis_list, np.array([[0, 0, 0]])),
     )
+
+    cartesian_basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), "cartesian")
+    cartesian_basis_list = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), ["cartesian"])
     assert np.allclose(
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type="cartesian"),
-        evaluate_basis(basis, np.array([[0, 0, 0]]), coord_type=["cartesian"]),
+        evaluate_basis(cartesian_basis, np.array([[0, 0, 0]])),
+        evaluate_basis(cartesian_basis_list, np.array([[0, 0, 0]])),
     )
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    spherical_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
+    spherical_basis_list = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8)
     assert np.allclose(
-        evaluate_basis(basis, np.array([[1, 1, 1]]), coord_type="spherical"),
-        evaluate_basis(basis, np.array([[1, 1, 1]]), coord_type=["spherical"] * 8),
+        evaluate_basis(spherical_basis, np.array([[1, 1, 1]])),
+        evaluate_basis(spherical_basis_list, np.array([[1, 1, 1]])),
     )
+
+    cartesian_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian")
+    cartesian_basis_list = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["cartesian"] * 8)
     assert np.allclose(
-        evaluate_basis(basis, np.array([[1, 1, 1]]), coord_type="cartesian"),
-        evaluate_basis(basis, np.array([[1, 1, 1]]), coord_type=["cartesian"] * 8),
+        evaluate_basis(cartesian_basis, np.array([[1, 1, 1]])),
+        evaluate_basis(cartesian_basis_list, np.array([[1, 1, 1]])),
     )
 
 
 def test_evaluate_basis_lincomb():
     """Test gbasis.evals.eval.evaluate_basis_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     evaluate_obj = Eval(basis)
     transform = np.random.rand(14, 18)
     assert np.allclose(
-        evaluate_obj.construct_array_lincomb(transform, "spherical", points=np.array([[1, 1, 1]])),
-        evaluate_basis(basis, np.array([[1, 1, 1]]), transform=transform, coord_type="spherical"),
+        evaluate_obj.construct_array_lincomb(transform, ["spherical"], points=np.array([[1, 1, 1]])),
+        evaluate_basis(basis, np.array([[1, 1, 1]]), transform=transform),
     )
 
 
@@ -134,7 +142,8 @@ def test_evaluate_basis_horton():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     points = np.array([[0, 0, 0], [0.8, 0, 0]])
     basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    cartesian_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, 'cartesian') for i in basis]
+    spherical_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, 'spherical') for i in basis]
 
     horton_eval_cart = np.load(find_datafile("data_horton_hhe_cart_eval.npy"))
     horton_eval_sph = np.load(find_datafile("data_horton_hhe_sph_eval.npy"))
@@ -143,8 +152,8 @@ def test_evaluate_basis_horton():
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
     grid_3d = np.vstack([grid_x.ravel(), grid_y.ravel(), grid_z.ravel()]).T
 
-    assert np.allclose(evaluate_basis(basis, grid_3d, coord_type="cartesian"), horton_eval_cart.T)
-    assert np.allclose(evaluate_basis(basis, grid_3d, coord_type="spherical"), horton_eval_sph.T)
+    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d), horton_eval_cart.T)
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d), horton_eval_sph.T)
 
 
 def test_evaluate_basis_pyscf():
diff --git a/tests/test_eval_deriv.py b/tests/test_eval_deriv.py
index e38c7a54..e00c07c1 100644
--- a/tests/test_eval_deriv.py
+++ b/tests/test_eval_deriv.py
@@ -131,14 +131,14 @@ def test_evaluate_deriv_construct_array_contraction():
 def test_evaluate_deriv_basis_cartesian():
     """Test gbasis.evals.eval.evaluate_deriv_basis_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian")
     evaluate_obj = EvalDeriv(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(
             points=np.array([[1, 1, 1]]), orders=np.array([0, 0, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), orders=np.array([0, 0, 0]), coord_type="cartesian"
+            basis, np.array([[1, 1, 1]]), orders=np.array([0, 0, 0])
         ),
     )
     assert np.allclose(
@@ -146,7 +146,7 @@ def test_evaluate_deriv_basis_cartesian():
             points=np.array([[1, 1, 1]]), orders=np.array([2, 1, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), orders=np.array([2, 1, 0]), coord_type="cartesian"
+            basis, np.array([[1, 1, 1]]), orders=np.array([2, 1, 0])
         ),
     )
 
@@ -154,14 +154,14 @@ def test_evaluate_deriv_basis_cartesian():
 def test_evaluate_deriv_basis_spherical():
     """Test gbasis.evals.eval.evaluate_deriv_basis_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
     evaluate_obj = EvalDeriv(basis)
     assert np.allclose(
         evaluate_obj.construct_array_spherical(
             points=np.array([[1, 1, 1]]), orders=np.array([0, 0, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), np.array([0, 0, 0]), coord_type="spherical"
+            basis, np.array([[1, 1, 1]]), np.array([0, 0, 0])
         ),
     )
     assert np.allclose(
@@ -169,7 +169,7 @@ def test_evaluate_deriv_basis_spherical():
             points=np.array([[1, 1, 1]]), orders=np.array([2, 1, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), np.array([2, 1, 0]), coord_type="spherical"
+            basis, np.array([[1, 1, 1]]), np.array([2, 1, 0])
         ),
     )
 
@@ -177,22 +177,24 @@ def test_evaluate_deriv_basis_spherical():
 def test_evaluate_deriv_basis_mix():
     """Test gbasis.evals.eval.evaluate_deriv_basis_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
-    evaluate_obj = EvalDeriv(basis)
+    cartesian_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["cartesian"] * 8)
+    spherical_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8)
+    evaluate_obj_cartesian = EvalDeriv(cartesian_basis)
+    evaluate_obj_spherical = EvalDeriv(spherical_basis)
     assert np.allclose(
-        evaluate_obj.construct_array_mix(
+        evaluate_obj_cartesian.construct_array_mix(
             ["cartesian"] * 8, points=np.array([[1, 1, 1]]), orders=np.array([0, 0, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), np.array([0, 0, 0]), coord_type=["cartesian"] * 8
+            cartesian_basis, np.array([[1, 1, 1]]), np.array([0, 0, 0])
         ),
     )
     assert np.allclose(
-        evaluate_obj.construct_array_mix(
+        evaluate_obj_spherical.construct_array_mix(
             ["spherical"] * 8, points=np.array([[1, 1, 1]]), orders=np.array([2, 1, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), np.array([2, 1, 0]), coord_type=["spherical"] * 8
+            spherical_basis, np.array([[1, 1, 1]]), np.array([2, 1, 0])
         ),
     )
 
@@ -200,27 +202,28 @@ def test_evaluate_deriv_basis_mix():
 def test_evaluate_deriv_basis_lincomb():
     """Test gbasis.evals.eval.evaluate_deriv_basis_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
-    evaluate_obj = EvalDeriv(basis)
+    cartesian_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian")
+    spherical_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
+    evaluate_obj_cartesian = EvalDeriv(cartesian_basis)
+    evaluate_obj_spherical = EvalDeriv(spherical_basis)
     cart_transform = np.random.rand(14, 19)
     sph_transform = np.random.rand(14, 18)
     assert np.allclose(
-        evaluate_obj.construct_array_lincomb(
-            cart_transform, "cartesian", points=np.array([[1, 1, 1]]), orders=np.array([0, 0, 0])
+        evaluate_obj_cartesian.construct_array_lincomb(
+            cart_transform, ["cartesian"], points=np.array([[1, 1, 1]]), orders=np.array([0, 0, 0])
         ),
         evaluate_deriv_basis(
-            basis,
+            cartesian_basis,
             np.array([[1, 1, 1]]),
             np.array([0, 0, 0]),
             cart_transform,
-            coord_type="cartesian",
         ),
     )
     assert np.allclose(
-        evaluate_obj.construct_array_lincomb(
-            sph_transform, "spherical", points=np.array([[1, 1, 1]]), orders=np.array([2, 1, 0])
+        evaluate_obj_spherical.construct_array_lincomb(
+            sph_transform, ["spherical"], points=np.array([[1, 1, 1]]), orders=np.array([2, 1, 0])
         ),
         evaluate_deriv_basis(
-            basis, np.array([[1, 1, 1]]), np.array([2, 1, 0]), sph_transform, coord_type="spherical"
+            spherical_basis, np.array([[1, 1, 1]]), np.array([2, 1, 0]), sph_transform
         ),
     )
diff --git a/tests/test_kinetic_energy.py b/tests/test_kinetic_energy.py
index 0abcf005..7bad067c 100644
--- a/tests/test_kinetic_energy.py
+++ b/tests/test_kinetic_energy.py
@@ -160,11 +160,11 @@ def test_kinetic_energy_construct_array_contraction():
 def test_kinetic_energy_integral_cartesian():
     """Test gbasis.integrals.kinetic_energy.kinetic_energy_integral_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     kinetic_energy_integral_obj = KineticEnergyIntegral(basis)
     assert np.allclose(
         kinetic_energy_integral_obj.construct_array_cartesian(),
-        kinetic_energy_integral(basis, coord_type="cartesian"),
+        kinetic_energy_integral(basis),
     )
 
 
@@ -172,11 +172,11 @@ def test_kinetic_energy_integral_spherical():
     """Test gbasis.integrals.kinetic_energy.kinetic_energy_integral_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     kinetic_energy_integral_obj = KineticEnergyIntegral(basis)
     assert np.allclose(
         kinetic_energy_integral_obj.construct_array_spherical(),
-        kinetic_energy_integral(basis, coord_type="spherical"),
+        kinetic_energy_integral(basis),
     )
 
 
@@ -184,23 +184,23 @@ def test_kinetic_energy_integral_mix():
     """Test gbasis.integrals.kinetic_energy.kinetic_energy_integral_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ['spherical'] * 8)
     kinetic_energy_integral_obj = KineticEnergyIntegral(basis)
     assert np.allclose(
         kinetic_energy_integral_obj.construct_array_mix(["spherical"] * 8),
-        kinetic_energy_integral(basis, coord_type=["spherical"] * 8),
+        kinetic_energy_integral(basis),
     )
 
 
 def test_kinetic_energy_integral_lincomb():
     """Test gbasis.integrals.kinetic_energy.kinetic_energy_integral_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     kinetic_energy_integral_obj = KineticEnergyIntegral(basis)
     transform = np.random.rand(14, 18)
     assert np.allclose(
-        kinetic_energy_integral_obj.construct_array_lincomb(transform, "spherical"),
-        kinetic_energy_integral(basis, transform, coord_type="spherical"),
+        kinetic_energy_integral_obj.construct_array_lincomb(transform, ["spherical"]),
+        kinetic_energy_integral(basis, transform),
     )
 
 
@@ -213,15 +213,15 @@ def test_kinetic_energy_integral_horton_anorcc_hhe():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     basis = make_contractions(
-        basis_dict, ["H", "He"], np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
+        basis_dict, ["H", "He"], np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]]), 'cartesian'
     )
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_kinetic_energy_integral = np.load(
         find_datafile("data_horton_hhe_cart_kinetic_energy_integral.npy")
     )
     assert np.allclose(
-        kinetic_energy_integral(basis, coord_type="cartesian"), horton_kinetic_energy_integral
+        kinetic_energy_integral(basis), horton_kinetic_energy_integral
     )
 
 
@@ -234,13 +234,13 @@ def test_kinetic_energy_integral_horton_anorcc_bec():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     basis = make_contractions(
-        basis_dict, ["Be", "C"], np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]])
+        basis_dict, ["Be", "C"], np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]]), 'cartesian'
     )
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_kinetic_energy_integral = np.load(
         find_datafile("data_horton_bec_cart_kinetic_energy_integral.npy")
     )
     assert np.allclose(
-        kinetic_energy_integral(basis, coord_type="cartesian"), horton_kinetic_energy_integral
+        kinetic_energy_integral(basis), horton_kinetic_energy_integral
     )
diff --git a/tests/test_moment.py b/tests/test_moment.py
index d57a1a66..2195f828 100644
--- a/tests/test_moment.py
+++ b/tests/test_moment.py
@@ -134,7 +134,7 @@ def test_moment_construct_array_contraction():
 def test_moment_cartesian():
     """Test gbasis.integrals.moment.moment_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     moment_obj = Moment(basis)
 
     assert np.allclose(
@@ -172,7 +172,6 @@ def test_moment_cartesian():
                     [0, 1, 1],
                 ]
             ),
-            coord_type="cartesian",
         ),
     )
 
@@ -181,7 +180,7 @@ def test_moment_spherical():
     """Test gbasis.integrals.moment.moment_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     moment_obj = Moment(basis)
     assert np.allclose(
         moment_obj.construct_array_spherical(
@@ -218,7 +217,6 @@ def test_moment_spherical():
                     [0, 1, 1],
                 ]
             ),
-            coord_type="spherical",
         ),
     )
 
@@ -227,7 +225,7 @@ def test_moment_mix():
     """Test gbasis.integrals.moment.moment_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ['spherical'] * 8)
     moment_obj = Moment(basis)
     assert np.allclose(
         moment_obj.construct_array_mix(
@@ -265,7 +263,6 @@ def test_moment_mix():
                     [0, 1, 1],
                 ]
             ),
-            coord_type=["spherical"] * 8,
         ),
     )
 
@@ -279,7 +276,7 @@ def test_moment_spherical_lincomb():
     assert np.allclose(
         moment_obj.construct_array_lincomb(
             transform,
-            "spherical",
+            ["spherical"],
             moment_coord=np.zeros(3),
             moment_orders=np.array(
                 [
diff --git a/tests/test_momentum.py b/tests/test_momentum.py
index b87fd3bd..5bd2c865 100644
--- a/tests/test_momentum.py
+++ b/tests/test_momentum.py
@@ -162,11 +162,11 @@ def test_momentum_construct_array_contraction():
 def test_momentum_integral_cartesian():
     """Test gbasis.integrals.momentum.momentum_integral_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     momentum_integral_obj = MomentumIntegral(basis)
     assert np.allclose(
         momentum_integral_obj.construct_array_cartesian(),
-        momentum_integral(basis, coord_type="cartesian"),
+        momentum_integral(basis),
     )
 
 
@@ -174,11 +174,11 @@ def test_momentum_integral_spherical():
     """Test gbasis.integrals.momentum.momentum_integral_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     momentum_integral_obj = MomentumIntegral(basis)
     assert np.allclose(
         momentum_integral_obj.construct_array_spherical(),
-        momentum_integral(basis, coord_type="spherical"),
+        momentum_integral(basis),
     )
 
 
@@ -186,21 +186,21 @@ def test_momentum_integral_mix():
     """Test gbasis.integrals.momentum.momentum_integral_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ['spherical'] * 8)
     momentum_integral_obj = MomentumIntegral(basis)
     assert np.allclose(
         momentum_integral_obj.construct_array_mix(["spherical"] * 8),
-        momentum_integral(basis, coord_type=["spherical"] * 8),
+        momentum_integral(basis),
     )
 
 
 def test_momentum_integral_lincomb():
     """Test gbasis.integrals.momentum.momentum_integral_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     momentum_integral_obj = MomentumIntegral(basis)
     transform = np.random.rand(14, 18)
     assert np.allclose(
-        momentum_integral_obj.construct_array_lincomb(transform, "spherical"),
-        momentum_integral(basis, transform=transform, coord_type="spherical"),
+        momentum_integral_obj.construct_array_lincomb(transform, ["spherical"]),
+        momentum_integral(basis, transform=transform),
     )
diff --git a/tests/test_nuclear_electron_attraction.py b/tests/test_nuclear_electron_attraction.py
index e306bdb0..2bbf00aa 100644
--- a/tests/test_nuclear_electron_attraction.py
+++ b/tests/test_nuclear_electron_attraction.py
@@ -15,13 +15,13 @@ def test_nuclear_electron_attraction_horton_anorcc_hhe():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], coords, 'cartesian')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_nucattract = np.load(find_datafile("data_horton_hhe_cart_nucattract.npy"))
     assert np.allclose(
         nuclear_electron_attraction_integral(
-            basis, coords, np.array([1, 2]), coord_type="cartesian"
+            basis, coords, np.array([1, 2])
         ),
         horton_nucattract,
     )
@@ -36,13 +36,13 @@ def test_nuclear_electron_attraction_horton_anorcc_bec():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     coords = np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["Be", "C"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["Be", "C"], coords, 'cartesian')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_nucattract = np.load(find_datafile("data_horton_bec_cart_nucattract.npy"))
     assert np.allclose(
         nuclear_electron_attraction_integral(
-            basis, coords, np.array([4, 6]), coord_type="cartesian"
+            basis, coords, np.array([4, 6])
         ),
         horton_nucattract,
     )
@@ -51,15 +51,15 @@ def test_nuclear_electron_attraction_horton_anorcc_bec():
 def test_nuclear_electron_attraction_cartesian():
     """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
 
     nuclear_coords = np.random.rand(5, 3)
     nuclear_charges = np.random.rand(5)
-    ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, coord_type="cartesian")
+    ref = point_charge_integral(basis, nuclear_coords, nuclear_charges)
     assert np.allclose(
         ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4],
         nuclear_electron_attraction_integral(
-            basis, nuclear_coords, nuclear_charges, coord_type="cartesian"
+            basis, nuclear_coords, nuclear_charges
         ),
     )
 
@@ -67,15 +67,15 @@ def test_nuclear_electron_attraction_cartesian():
 def test_nuclear_electron_attraction_spherical():
     """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
 
     nuclear_coords = np.random.rand(5, 3)
     nuclear_charges = np.random.rand(5)
-    ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, coord_type="spherical")
+    ref = point_charge_integral(basis, nuclear_coords, nuclear_charges)
     assert np.allclose(
         ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4],
         nuclear_electron_attraction_integral(
-            basis, nuclear_coords, nuclear_charges, coord_type="spherical"
+            basis, nuclear_coords, nuclear_charges
         ),
     )
 
@@ -83,17 +83,17 @@ def test_nuclear_electron_attraction_spherical():
 def test_nuclear_electron_attraction_mix():
     """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ['spherical'] * 8)
 
     nuclear_coords = np.random.rand(5, 3)
     nuclear_charges = np.random.rand(5)
     ref = point_charge_integral(
-        basis, nuclear_coords, nuclear_charges, coord_type=["spherical"] * 8
+        basis, nuclear_coords, nuclear_charges
     )
     assert np.allclose(
         ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4],
         nuclear_electron_attraction_integral(
-            basis, nuclear_coords, nuclear_charges, coord_type=["spherical"] * 8
+            basis, nuclear_coords, nuclear_charges
         ),
     )
 
@@ -101,7 +101,7 @@ def test_nuclear_electron_attraction_mix():
 def test_nuclear_electron_attraction_lincomb():
     """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]),)
 
     nuclear_coords = np.random.rand(5, 3)
     nuclear_charges = np.random.rand(5)
diff --git a/tests/test_overlap.py b/tests/test_overlap.py
index 59ca1346..8a51fcf9 100644
--- a/tests/test_overlap.py
+++ b/tests/test_overlap.py
@@ -75,10 +75,10 @@ def test_overlap_construct_array_contraction():
 def test_overlap_cartesian():
     """Test gbasis.integrals.overlap.overlap_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian")
     overlap_obj = Overlap(basis)
     assert np.allclose(
-        overlap_obj.construct_array_cartesian(), overlap_integral(basis, coord_type="cartesian")
+        overlap_obj.construct_array_cartesian(), overlap_integral(basis)
     )
 
 
@@ -86,10 +86,10 @@ def test_overlap_spherical():
     """Test gbasis.integrals.overlap.overlap_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
     overlap_obj = Overlap(basis)
     assert np.allclose(
-        overlap_obj.construct_array_spherical(), overlap_integral(basis, coord_type="spherical")
+        overlap_obj.construct_array_spherical(), overlap_integral(basis)
     )
 
 
@@ -97,23 +97,23 @@ def test_overlap_mix():
     """Test gbasis.integrals.overlap.overlap_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8)
     overlap_obj = Overlap(basis)
     assert np.allclose(
         overlap_obj.construct_array_mix(["spherical"] * 8),
-        overlap_integral(basis, coord_type=["spherical"] * 8),
+        overlap_integral(basis),
     )
 
 
 def test_overlap_lincomb():
     """Test gbasis.integrals.overlap.overlap_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
     overlap_obj = Overlap(basis)
     transform = np.random.rand(14, 18)
     assert np.allclose(
-        overlap_obj.construct_array_lincomb(transform, "spherical"),
-        overlap_integral(basis, transform=transform, coord_type="spherical"),
+        overlap_obj.construct_array_lincomb(transform, ["spherical"]),
+        overlap_integral(basis, transform=transform),
     )
 
 
@@ -170,7 +170,7 @@ def test_overlap_cartesian_norm_sto6g():
     """
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     overlap_obj = Overlap(basis)
     assert np.allclose(np.diag(overlap_obj.construct_array_cartesian()), 1)
 
@@ -184,12 +184,12 @@ def test_overlap_horton_anorcc_hhe():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     basis = make_contractions(
-        basis_dict, ["H", "He"], np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
+        basis_dict, ["H", "He"], np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]]), "cartesian"
     )
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_overlap = np.load(find_datafile("data_horton_hhe_cart_overlap.npy"))
-    assert np.allclose(overlap_integral(basis, coord_type="cartesian"), horton_overlap)
+    assert np.allclose(overlap_integral(basis), horton_overlap)
 
 
 def test_overlap_horton_anorcc_bec():
@@ -201,9 +201,10 @@ def test_overlap_horton_anorcc_bec():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     basis = make_contractions(
-        basis_dict, ["Be", "C"], np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]])
+        basis_dict, ["Be", "C"], np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]]), "cartesian"
     )
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_overlap = np.load(find_datafile("data_horton_bec_cart_overlap.npy"))
-    assert np.allclose(overlap_integral(basis, coord_type="cartesian"), horton_overlap)
+    assert np.allclose(overlap_integral(basis), horton_overlap)
diff --git a/tests/test_overlap_asymm.py b/tests/test_overlap_asymm.py
index d2cd20f5..e6f43529 100644
--- a/tests/test_overlap_asymm.py
+++ b/tests/test_overlap_asymm.py
@@ -15,14 +15,14 @@ def test_overlap_integral_asymmetric_horton_anorcc_hhe():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     basis = make_contractions(
-        basis_dict, ["H", "He"], np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
+        basis_dict, ["H", "He"], np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]]), "cartesian"
     )
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_overlap_integral_asymmetric = np.load(find_datafile("data_horton_hhe_cart_overlap.npy"))
     assert np.allclose(
         overlap_integral_asymmetric(
-            basis, basis, coord_type_two="cartesian", coord_type_one="cartesian"
+            basis, basis
         ),
         horton_overlap_integral_asymmetric,
     )
@@ -37,14 +37,14 @@ def test_overlap_integral_asymmetric_horton_anorcc_bec():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     basis = make_contractions(
-        basis_dict, ["Be", "C"], np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]])
+        basis_dict, ["Be", "C"], np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]]), "cartesian"
     )
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_overlap_integral_asymmetric = np.load(find_datafile("data_horton_bec_cart_overlap.npy"))
     assert np.allclose(
         overlap_integral_asymmetric(
-            basis, basis, coord_type_two="cartesian", coord_type_one="cartesian"
+            basis, basis
         ),
         horton_overlap_integral_asymmetric,
     )
@@ -53,39 +53,38 @@ def test_overlap_integral_asymmetric_horton_anorcc_bec():
 def test_overlap_integral_asymmetric_compare():
     """Test overlap_asymm.overlap_integral_asymmetric against overlap.overlap_integral."""
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
-
-    basis = make_contractions(basis_dict, ["Kr", "Kr"], np.array([[0, 0, 0], [1.0, 0, 0]]))
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    cartesian_basis = make_contractions(basis_dict, ["Kr", "Kr"], np.array([[0, 0, 0], [1.0, 0, 0]]), "cartesian")
+    spherical_basis = make_contractions(basis_dict, ["Kr", "Kr"], np.array([[0, 0, 0], [1.0, 0, 0]]), "spherical")
+    mixed_basis = make_contractions(basis_dict, ["Kr", "Kr"], np.array([[0, 0, 0], [1.0, 0, 0]]), ['spherical'] * 9 + ['cartesian'])
+    cartesian_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in cartesian_basis]
+    spherical_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in spherical_basis]
+    mixed_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in mixed_basis]
 
     assert np.allclose(
-        overlap_integral(basis, coord_type="cartesian"),
+        overlap_integral(cartesian_basis),
         overlap_integral_asymmetric(
-            basis, basis, coord_type_one="cartesian", coord_type_two="cartesian"
+            cartesian_basis, cartesian_basis
         ),
     )
     assert np.allclose(
-        overlap_integral(basis, coord_type="spherical"),
+        overlap_integral(spherical_basis),
         overlap_integral_asymmetric(
-            basis, basis, coord_type_one="spherical", coord_type_two="spherical"
+            spherical_basis, spherical_basis
         ),
     )
     assert np.allclose(
-        overlap_integral(basis, transform=np.identity(218), coord_type="spherical"),
+        overlap_integral(spherical_basis, transform=np.identity(218)),
         overlap_integral_asymmetric(
-            basis,
-            basis,
+            spherical_basis,
+            spherical_basis,
             transform_one=np.identity(218),
             transform_two=np.identity(218),
-            coord_type_one="spherical",
-            coord_type_two="spherical",
         ),
     )
     assert np.allclose(
-        overlap_integral(basis, coord_type=["spherical"] * 9 + ["cartesian"]),
+        overlap_integral(mixed_basis),
         overlap_integral_asymmetric(
-            basis,
-            basis,
-            coord_type_one=["spherical"] * 9 + ["cartesian"],
-            coord_type_two=["spherical"] * 9 + ["cartesian"],
+            mixed_basis,
+            mixed_basis,
         ),
     )
diff --git a/tests/test_point_charge.py b/tests/test_point_charge.py
index f79da5a8..3add66ce 100644
--- a/tests/test_point_charge.py
+++ b/tests/test_point_charge.py
@@ -119,7 +119,7 @@ def test_construct_array_contraction():
 def test_point_charge_cartesian():
     """Test gbasis.integrals.point_charge.point_charge_cartesian."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     point_charge_obj = PointChargeIntegral(basis)
 
     points_coords = np.random.rand(5, 3)
@@ -129,7 +129,7 @@ def test_point_charge_cartesian():
             points_coords=points_coords, points_charge=points_charge
         ),
         point_charge_integral(
-            basis, points_coords=points_coords, points_charge=points_charge, coord_type="cartesian"
+            basis, points_coords=points_coords, points_charge=points_charge
         ),
     )
 
@@ -137,7 +137,7 @@ def test_point_charge_cartesian():
 def test_point_charge_spherical():
     """Test gbasis.integrals.point_charge.point_charge_spherical."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     point_charge_obj = PointChargeIntegral(basis)
 
     points_coords = np.random.rand(5, 3)
@@ -147,7 +147,7 @@ def test_point_charge_spherical():
             points_coords=points_coords, points_charge=points_charge
         ),
         point_charge_integral(
-            basis, points_coords=points_coords, points_charge=points_charge, coord_type="spherical"
+            basis, points_coords=points_coords, points_charge=points_charge
         ),
     )
 
@@ -155,7 +155,7 @@ def test_point_charge_spherical():
 def test_point_charge_mix():
     """Test gbasis.integrals.point_charge.point_charge_mix."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8)
     point_charge_obj = PointChargeIntegral(basis)
 
     points_coords = np.random.rand(5, 3)
@@ -164,14 +164,14 @@ def test_point_charge_mix():
         point_charge_obj.construct_array_mix(
             ["spherical"] * 8, points_coords=points_coords, points_charge=points_charge
         ),
-        point_charge_integral(basis, points_coords, points_charge, coord_type=["spherical"] * 8),
+        point_charge_integral(basis, points_coords, points_charge)
     )
 
 
 def test_point_charge_lincomb():
     """Test gbasis.integrals.point_charge.point_charge_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     point_charge_obj = PointChargeIntegral(basis)
 
     points_coords = np.random.rand(5, 3)
@@ -179,7 +179,7 @@ def test_point_charge_lincomb():
     transform = np.random.rand(14, 18)
     assert np.allclose(
         point_charge_obj.construct_array_lincomb(
-            transform, "spherical", points_coords=points_coords, points_charge=points_charge
+            transform, ["spherical"], points_coords=points_coords, points_charge=points_charge
         ),
         point_charge_integral(
             basis, points_coords=points_coords, points_charge=points_charge, transform=transform
diff --git a/tests/test_wrappers.py b/tests/test_wrappers.py
index bb5aa30c..70c6e077 100644
--- a/tests/test_wrappers.py
+++ b/tests/test_wrappers.py
@@ -14,7 +14,8 @@ def test_from_iodata():
 
     mol = load_one(find_datafile("data_iodata_water_sto3g_hf_g03.fchk"))
 
-    basis, coord_types = from_iodata(mol)
+    basis = from_iodata(mol)
+    coord_types = [type for type in [shell.coord_type for shell in basis]]
 
     assert coord_types == ["cartesian"] * 5
     assert all(isinstance(i, GeneralizedContractionShell) for i in basis)
@@ -77,28 +78,30 @@ def test_from_iodata():
     # NOTE: you shouldn't actually change the magnetic quantum number that is not compatible with
     # the angular momentum, but we do so here to check that user input is accepted
     mol.obasis.conventions[(0, "p")] = ["c1"]
-    basis, coord_types = from_iodata(mol)
+    basis = from_iodata(mol)
+    coord_types = [type for type in [shell.coord_type for shell in basis]]
     basis[2].angmom = 0
     assert coord_types == ["cartesian"] * 5
     assert basis[2].angmom_components_sph == ("c1",)
     assert np.allclose(basis[2].norm_cont, 1.0)
 
     mol.obasis.conventions[(1, "p")] = ["c1", "c0", "s1"]
-    basis, coord_types = from_iodata(mol)
+    basis = from_iodata(mol)
+    coord_types = [type for type in [shell.coord_type for shell in basis]]
     basis[2].angmom = 1
     assert coord_types == ["cartesian"] * 5
     assert basis[2].angmom_components_sph == ("c1", "c0", "s1")
     assert np.allclose(basis[2].norm_cont, 1.0)
 
     mol.obasis.conventions[(1, "c")] = ["z", "y", "x"]
-    basis, coord_types = from_iodata(mol)
+    basis = from_iodata(mol)
     basis[2].angmom = 1
     assert np.allclose(np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]]), basis[2].angmom_components_cart)
 
     # Test Cartesian convention generation for missing angmom
     # Needed for cases when only spherical basis is used in basis set
     del mol.obasis.conventions[(1, "c")]
-    basis, coord_types = from_iodata(mol)
+    basis = from_iodata(mol)
     basis[2].angmom = 1
     assert np.allclose(np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), basis[2].angmom_components_cart)
 
@@ -108,7 +111,7 @@ def test_from_iodata():
 
     with pytest.raises(ValueError):
         mol.obasis = mol.obasis._replace(primitive_normalization="L1")
-        basis, coord_types = from_iodata(mol)
+        basis = from_iodata(mol)
 
 
 def test_from_pyscf():

From 5a52509a1ec4cf32bfe79573439aadd5489825f0 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Sun, 24 Dec 2023 23:46:28 +0300
Subject: [PATCH 07/17] Fix TypeError in test_evaluate_basis_pyscf

Fix the error below:
evaluate_basis() got an unexpected keyword argument 'coord_type'
---
 tests/test_eval.py | 48 ++++++++++++----------------------------------
 1 file changed, 12 insertions(+), 36 deletions(-)

diff --git a/tests/test_eval.py b/tests/test_eval.py
index 1765c5b4..70e6b4ac 100644
--- a/tests/test_eval.py
+++ b/tests/test_eval.py
@@ -175,45 +175,21 @@ def test_evaluate_basis_pyscf():
     pyscf_eval_cart = gto.eval_gto(mol, "GTOval_cart", grid_3d)
 
     # s orbitals
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[:6], pyscf_eval_cart.T[:6]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[46:53], pyscf_eval_cart.T[46:53]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[:6], pyscf_eval_sph.T[:6]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[40:47], pyscf_eval_sph.T[40:47]
-    )
+    assert np.allclose(evaluate_basis(basis, grid_3d)[:6], pyscf_eval_cart.T[:6])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[46:53], pyscf_eval_cart.T[46:53])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[:6], pyscf_eval_sph.T[:6])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[40:47], pyscf_eval_sph.T[40:47])
     # p orbitals
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[6:18], pyscf_eval_cart.T[6:18]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[53:65], pyscf_eval_cart.T[53:65]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[6:18], pyscf_eval_sph.T[6:18]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[47:59], pyscf_eval_sph.T[47:59]
-    )
+    assert np.allclose(evaluate_basis(basis, grid_3d)[6:18], pyscf_eval_cart.T[6:18])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[53:65], pyscf_eval_cart.T[53:65])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[6:18], pyscf_eval_sph.T[6:18])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[47:59], pyscf_eval_sph.T[47:59])
     # d orbitals are off by some constant for the cartesian case
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[18:33], pyscf_eval_sph.T[18:33]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[59:74], pyscf_eval_sph.T[59:74]
-    )
+    assert np.allclose(evaluate_basis(basis, grid_3d)[18:33], pyscf_eval_sph.T[18:33])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[59:74], pyscf_eval_sph.T[59:74])
     # f orbitals are off by some constant for the cartesian case
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[33:40], pyscf_eval_sph.T[33:40]
-    )
-    assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="spherical")[74:88], pyscf_eval_sph.T[74:88]
-    )
+    assert np.allclose(evaluate_basis(basis, grid_3d)[33:40], pyscf_eval_sph.T[33:40])
+    assert np.allclose(evaluate_basis(basis, grid_3d)[74:88], pyscf_eval_sph.T[74:88])
 
 
 @pytest.mark.xfail

From 6109e8d850a2e41e025cadaf49d7ab06ce3349b0 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Sun, 24 Dec 2023 23:05:50 +0300
Subject: [PATCH 08/17] Add options for coord_type to docstring

---
 gbasis/contractions.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/gbasis/contractions.py b/gbasis/contractions.py
index 06887495..42a34564 100644
--- a/gbasis/contractions.py
+++ b/gbasis/contractions.py
@@ -116,7 +116,7 @@ class GeneralizedContractionShell:
 
     """
 
-    def __init__(self, angmom, coord, coeffs, exps, coord_type='p'):
+    def __init__(self, angmom, coord, coeffs, exps, coord_type="p"):
         r"""Initialize a GeneralizedContractionShell instance.
 
         Parameters
@@ -137,9 +137,9 @@ def __init__(self, angmom, coord, coeffs, exps, coord_type='p'):
             dimension.
         exps : np.ndarray(K,)
             Exponents of the primitives, :math:`\{\alpha_i\}_{i=1}^K`.
-        coord_type : str
-            Coordinate type of the contraction.
-            Default is 'p' (spherical).
+        coord_type : str, optional
+            Coordinate type of the contraction. Options include "cartesian" or "c" and
+            "spherical" or "p".
 
         """
         self.angmom = angmom

From 155c1a96c6b4cd30cdd9f29358a8a1aa28551651 Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Tue, 9 Jan 2024 15:22:59 -0500
Subject: [PATCH 09/17] Fix failed PySCF test

I took a similar approach to `make_contractions` for creating
a `basis` with PySCF - now the user specifies the coordinate type of each
shell using either a string or a list of strings
---
 gbasis/wrappers.py | 27 +++++++++++++++++++++++++--
 tests/test_eval.py | 27 ++++++++++++++-------------
 2 files changed, 39 insertions(+), 15 deletions(-)

diff --git a/gbasis/wrappers.py b/gbasis/wrappers.py
index b75cdffc..016db858 100644
--- a/gbasis/wrappers.py
+++ b/gbasis/wrappers.py
@@ -138,13 +138,19 @@ def angmom_components_sph(self):
     return basis
 
 
-def from_pyscf(mol):
+def from_pyscf(mol, coord_types="spherical"):
     """Return basis set stored within the `Mole` instance in `pyscf`.
 
     Parameters
     ----------
     mol : pyscf.gto.mole.Mole
         `Mole` object in `pyscf`.
+    coord_types : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
+        Types of the coordinate system for the contractions.
+        If "cartesian", then all of the contractions are treated as Cartesian contractions.
+        If "spherical", then all of the contractions are treated as spherical contractions.
+        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
+        coordinate type of each `PyscfShell` instance.
 
     Returns
     -------
@@ -183,6 +189,7 @@ def angmom_components_sph(self):
             return super().angmom_components_sph
 
     basis = []
+    len_coord_types = len(coord_types)
     for atom, coord in mol._atom:
         basis_info = mol._basis[atom]
 
@@ -194,6 +201,22 @@ def angmom_components_sph(self):
 
             coeffs = exps_coeffs[:, 1:]
 
-            basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps)))
+            if type(coord_types) == str:
+            # if coord_types given as a single string, assign the specified type to all contractions for all atoms
+                if coord_types == "spherical":
+                    basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps), 'p'))
+                elif coord_types == "cartesian":
+                    basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps), 'c'))
+                else:
+                    raise ValueError("If coord_types is a string, it must be either 'spherical' or 'cartesian'.")
+            elif type(coord_types) == list:
+                # if coord_types given as a list, assign the specified type to each atom's contractions individually
+                if len_coord_types == mol.nbas:
+                    basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps), coord_types.pop(0)))
+                else:
+                    raise ValueError(
+                        "If coord_types is a list, it must be the same length as the total number of contractions.")
+            else:
+                raise TypeError("coord_types must be a string or list of strings.")
 
     return tuple(basis)
diff --git a/tests/test_eval.py b/tests/test_eval.py
index 70e6b4ac..a0f28711 100644
--- a/tests/test_eval.py
+++ b/tests/test_eval.py
@@ -165,7 +165,8 @@ def test_evaluate_basis_pyscf():
 
     mol = gto.Mole()
     mol.build(atom="H 0 0 0; He 0.8 0 0", basis="ano-rcc", spin=1)
-    basis = from_pyscf(mol)
+    cartesian_basis = from_pyscf(mol, coord_types="cartesian")
+    spherical_basis = from_pyscf(mol, coord_types="spherical")
 
     grid_1d = np.linspace(-2, 2, num=5)
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
@@ -175,21 +176,21 @@ def test_evaluate_basis_pyscf():
     pyscf_eval_cart = gto.eval_gto(mol, "GTOval_cart", grid_3d)
 
     # s orbitals
-    assert np.allclose(evaluate_basis(basis, grid_3d)[:6], pyscf_eval_cart.T[:6])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[46:53], pyscf_eval_cart.T[46:53])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[:6], pyscf_eval_sph.T[:6])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[40:47], pyscf_eval_sph.T[40:47])
+    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[:6], pyscf_eval_cart.T[:6])
+    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[46:53], pyscf_eval_cart.T[46:53])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[:6], pyscf_eval_sph.T[:6])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[40:47], pyscf_eval_sph.T[40:47])
     # p orbitals
-    assert np.allclose(evaluate_basis(basis, grid_3d)[6:18], pyscf_eval_cart.T[6:18])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[53:65], pyscf_eval_cart.T[53:65])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[6:18], pyscf_eval_sph.T[6:18])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[47:59], pyscf_eval_sph.T[47:59])
+    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[6:18], pyscf_eval_cart.T[6:18])
+    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[53:65], pyscf_eval_cart.T[53:65])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[6:18], pyscf_eval_sph.T[6:18])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[47:59], pyscf_eval_sph.T[47:59])
     # d orbitals are off by some constant for the cartesian case
-    assert np.allclose(evaluate_basis(basis, grid_3d)[18:33], pyscf_eval_sph.T[18:33])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[59:74], pyscf_eval_sph.T[59:74])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[18:33], pyscf_eval_sph.T[18:33])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[59:74], pyscf_eval_sph.T[59:74])
     # f orbitals are off by some constant for the cartesian case
-    assert np.allclose(evaluate_basis(basis, grid_3d)[33:40], pyscf_eval_sph.T[33:40])
-    assert np.allclose(evaluate_basis(basis, grid_3d)[74:88], pyscf_eval_sph.T[74:88])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[33:40], pyscf_eval_sph.T[33:40])
+    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[74:88], pyscf_eval_sph.T[74:88])
 
 
 @pytest.mark.xfail

From d6a35db836b3e42dcb9f672de8ac31607d1f5c1e Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Tue, 9 Jan 2024 17:04:48 -0500
Subject: [PATCH 10/17] Stop using reserved type keyword improperly

---
 gbasis/evals/electrostatic_potential.py | 6 +++---
 gbasis/evals/eval.py                    | 6 +++---
 gbasis/evals/eval_deriv.py              | 6 +++---
 gbasis/integrals/angular_momentum.py    | 6 +++---
 gbasis/integrals/electron_repulsion.py  | 6 +++---
 gbasis/integrals/kinetic_energy.py      | 6 +++---
 gbasis/integrals/moment.py              | 6 +++---
 gbasis/integrals/momentum.py            | 6 +++---
 gbasis/integrals/overlap.py             | 6 +++---
 gbasis/integrals/overlap_asymm.py       | 4 ++--
 gbasis/integrals/point_charge.py        | 6 +++---
 11 files changed, 32 insertions(+), 32 deletions(-)

diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py
index 10b72b0c..d90dd3ca 100644
--- a/gbasis/evals/electrostatic_potential.py
+++ b/gbasis/evals/electrostatic_potential.py
@@ -93,15 +93,15 @@ def electrostatic_potential(
     if threshold_dist < 0:
         raise ValueError("`threshold_dist` must be greater than or equal to zero.")
 
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         if sum(cont.num_cart * cont.num_seg_cont for cont in basis) != one_density_matrix.shape[0]:
             raise ValueError(
                 "`one_density_matrix` does not have number of rows/columns that is equal to the "
                 "total number of Cartesian contractions (atomic orbitals)."
             )
-    elif all(type == "spherical" for type in coord_type):
+    elif all(ct == "spherical" for ct in coord_type):
         if sum(cont.num_sph * cont.num_seg_cont for cont in basis) != one_density_matrix.shape[0]:
             raise ValueError(
                 "`one_density_matrix` does not have number of rows/columns that is equal to the "
diff --git a/gbasis/evals/eval.py b/gbasis/evals/eval.py
index 54107012..bd14b8d6 100644
--- a/gbasis/evals/eval.py
+++ b/gbasis/evals/eval.py
@@ -140,12 +140,12 @@ def evaluate_basis(basis, points, transform=None):
         `N` is the number of coordinates at which the contractions are evaluated.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return Eval(basis).construct_array_lincomb(transform, coord_type, points=points)
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return Eval(basis).construct_array_cartesian(points=points)
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return Eval(basis).construct_array_spherical(points=points)
     return Eval(basis).construct_array_mix(coord_type, points=points)
diff --git a/gbasis/evals/eval_deriv.py b/gbasis/evals/eval_deriv.py
index a9cb7e6d..fd17343e 100644
--- a/gbasis/evals/eval_deriv.py
+++ b/gbasis/evals/eval_deriv.py
@@ -160,14 +160,14 @@ def evaluate_deriv_basis(basis, points, orders, transform=None):
         `N` is the number of coordinates at which the contractions are evaluated.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return EvalDeriv(basis).construct_array_lincomb(
             transform, coord_type, points=points, orders=orders
         )
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return EvalDeriv(basis).construct_array_cartesian(points=points, orders=orders)
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return EvalDeriv(basis).construct_array_spherical(points=points, orders=orders)
     return EvalDeriv(basis).construct_array_mix(coord_type, points=points, orders=orders)
diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py
index daefdd8a..084827d3 100644
--- a/gbasis/integrals/angular_momentum.py
+++ b/gbasis/integrals/angular_momentum.py
@@ -178,12 +178,12 @@ def angular_momentum_integral(basis, transform=None):
         Dimension 2 corresponds to the direction of the angular momentum :math:`(x, y, z)`.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return AngularMomentumIntegral(basis).construct_array_lincomb(transform, coord_type)
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return AngularMomentumIntegral(basis).construct_array_cartesian()
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return AngularMomentumIntegral(basis).construct_array_spherical()
     return AngularMomentumIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/electron_repulsion.py b/gbasis/integrals/electron_repulsion.py
index d3098de3..cbc575d3 100644
--- a/gbasis/integrals/electron_repulsion.py
+++ b/gbasis/integrals/electron_repulsion.py
@@ -238,13 +238,13 @@ def electron_repulsion_integral(
     if notation not in ["physicist", "chemist"]:
         raise ValueError("`notation` must be one of 'physicist' or 'chemist'")
 
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         array = ElectronRepulsionIntegral(basis).construct_array_lincomb(transform, coord_type)
-    elif all(type == "cartesian" for type in coord_type):
+    elif all(ct == "cartesian" for ct in coord_type):
         array = ElectronRepulsionIntegral(basis).construct_array_cartesian()
-    elif all(type == "spherical" for type in coord_type):
+    elif all(ct == "spherical" for ct in coord_type):
         array = ElectronRepulsionIntegral(basis).construct_array_spherical()
     else:
         array = ElectronRepulsionIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py
index 8ee4b8a3..d196bf02 100644
--- a/gbasis/integrals/kinetic_energy.py
+++ b/gbasis/integrals/kinetic_energy.py
@@ -142,12 +142,12 @@ def kinetic_energy_integral(basis, transform=None):
         `K_orbs` is the total number of basis functions in the basis set.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return KineticEnergyIntegral(basis).construct_array_lincomb(transform, coord_type)
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return KineticEnergyIntegral(basis).construct_array_cartesian()
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return KineticEnergyIntegral(basis).construct_array_spherical()
     return KineticEnergyIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py
index 6f1291fb..44ec36ef 100644
--- a/gbasis/integrals/moment.py
+++ b/gbasis/integrals/moment.py
@@ -194,17 +194,17 @@ def moment_integral(basis, moment_coord, moment_orders, transform=None):
     order.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return Moment(basis).construct_array_lincomb(
             transform, coord_type, moment_coord=moment_coord, moment_orders=moment_orders
         )
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return Moment(basis).construct_array_cartesian(
             moment_coord=moment_coord, moment_orders=moment_orders
         )
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return Moment(basis).construct_array_spherical(
             moment_coord=moment_coord, moment_orders=moment_orders
         )
diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py
index f9755c80..053d2fc6 100644
--- a/gbasis/integrals/momentum.py
+++ b/gbasis/integrals/momentum.py
@@ -133,12 +133,12 @@ def momentum_integral(basis, transform=None):
         number of basis functions in the basis set.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return MomentumIntegral(basis).construct_array_lincomb(transform, coord_type)
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return MomentumIntegral(basis).construct_array_cartesian()
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return MomentumIntegral(basis).construct_array_spherical()
     return MomentumIntegral(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py
index 10eae345..b9515bdb 100644
--- a/gbasis/integrals/overlap.py
+++ b/gbasis/integrals/overlap.py
@@ -129,12 +129,12 @@ def overlap_integral(basis, transform=None):
         number of basis functions in the basis set.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return Overlap(basis).construct_array_lincomb(transform, coord_type)
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return Overlap(basis).construct_array_cartesian()
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return Overlap(basis).construct_array_spherical()
     return Overlap(basis).construct_array_mix(coord_type)
diff --git a/gbasis/integrals/overlap_asymm.py b/gbasis/integrals/overlap_asymm.py
index 4f0858ee..38a78d4a 100644
--- a/gbasis/integrals/overlap_asymm.py
+++ b/gbasis/integrals/overlap_asymm.py
@@ -98,8 +98,8 @@ def overlap_integral_asymmetric(
         `K_orbs_2` is the number of basis functions in the `basis_two`.
 
     """
-    coord_type_one = [type for type in [shell.coord_type for shell in basis_one]]
-    coord_type_two = [type for type in [shell.coord_type for shell in basis_two]]
+    coord_type_one = [ct for ct in [shell.coord_type for shell in basis_one]]
+    coord_type_two = [ct for ct in [shell.coord_type for shell in basis_two]]
 
     return OverlapAsymmetric(basis_one, basis_two).construct_array_lincomb(
         transform_one, transform_two, coord_type_one, coord_type_two
diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py
index ab3d1857..55ef38af 100644
--- a/gbasis/integrals/point_charge.py
+++ b/gbasis/integrals/point_charge.py
@@ -297,17 +297,17 @@ def point_charge_integral(
         `N` is the number of coordinates at which the contractions are evaluated.
 
     """
-    coord_type = [type for type in [shell.coord_type for shell in basis]]
+    coord_type = [ct for ct in [shell.coord_type for shell in basis]]
 
     if transform is not None:
         return PointChargeIntegral(basis).construct_array_lincomb(
             transform, coord_type, points_coords=points_coords, points_charge=points_charge
         )
-    if all(type == "cartesian" for type in coord_type):
+    if all(ct == "cartesian" for ct in coord_type):
         return PointChargeIntegral(basis).construct_array_cartesian(
             points_coords=points_coords, points_charge=points_charge
         )
-    if all(type == "spherical" for type in coord_type):
+    if all(ct == "spherical" for ct in coord_type):
         return PointChargeIntegral(basis).construct_array_spherical(
             points_coords=points_coords, points_charge=points_charge
         )

From 729d83e92d19e87b545e5aff64b575639ea449b2 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Thu, 11 Jan 2024 01:31:04 +0330
Subject: [PATCH 11/17] Use pyscf.gto.Mole.cart to get coord_type

This makes it possible to simplify from_pyscf function to remove
coord_types argument.
---
 gbasis/wrappers.py | 12 +++++-------
 tests/test_eval.py | 11 ++++++++---
 2 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/gbasis/wrappers.py b/gbasis/wrappers.py
index 016db858..50c776be 100644
--- a/gbasis/wrappers.py
+++ b/gbasis/wrappers.py
@@ -138,19 +138,13 @@ def angmom_components_sph(self):
     return basis
 
 
-def from_pyscf(mol, coord_types="spherical"):
+def from_pyscf(mol):
     """Return basis set stored within the `Mole` instance in `pyscf`.
 
     Parameters
     ----------
     mol : pyscf.gto.mole.Mole
         `Mole` object in `pyscf`.
-    coord_types : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
-        Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
-        coordinate type of each `PyscfShell` instance.
 
     Returns
     -------
@@ -174,6 +168,10 @@ def from_pyscf(mol, coord_types="spherical"):
     if not (mol.__class__.__name__ == "Mole" and hasattr(mol, "_basis")):
         raise ValueError("`mol` must be a `pyscf.gto.mole.Mole` instance.")
 
+    # assign the coordinate types (which can be either Cartesian or Spherical)
+    # it seems like pyscf does not support mixed "cartesian" and "spherical" basis.
+    coord_types = "cartesian" if mol.cart else "spherical"
+
     class PyscfShell(GeneralizedContractionShell):
         """Shell object that is compatible with gbasis' shell object.
 
diff --git a/tests/test_eval.py b/tests/test_eval.py
index a0f28711..954477d1 100644
--- a/tests/test_eval.py
+++ b/tests/test_eval.py
@@ -163,10 +163,15 @@ def test_evaluate_basis_pyscf():
     from pyscf import gto
     from gbasis.wrappers import from_pyscf
 
-    mol = gto.Mole()
+    # build Cartesian basis (default is cart=False)
+    mol = gto.Mole(cart=True)
+    mol.build(atom="H 0 0 0; He 0.8 0 0", basis="ano-rcc", spin=1)
+    cartesian_basis = from_pyscf(mol)
+
+    # build Spherical basis (default is cart=False)
+    mol = gto.Mole(cart=False)
     mol.build(atom="H 0 0 0; He 0.8 0 0", basis="ano-rcc", spin=1)
-    cartesian_basis = from_pyscf(mol, coord_types="cartesian")
-    spherical_basis = from_pyscf(mol, coord_types="spherical")
+    spherical_basis = from_pyscf(mol)
 
     grid_1d = np.linspace(-2, 2, num=5)
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)

From db0bc54091f23038b0f30e0b00bf8c6d76c60447 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Thu, 11 Jan 2024 01:40:33 +0330
Subject: [PATCH 12/17] Replace variable name type which a built-in function
 name

---
 gbasis/base_four_symm.py | 4 ++--
 gbasis/base_one.py       | 4 ++--
 gbasis/base_two_asymm.py | 8 ++++++--
 gbasis/base_two_symm.py  | 4 ++--
 4 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/gbasis/base_four_symm.py b/gbasis/base_four_symm.py
index 03f8fa1c..6f69c336 100644
--- a/gbasis/base_four_symm.py
+++ b/gbasis/base_four_symm.py
@@ -595,9 +595,9 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
             If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if all(type == "cartesian" for type in coord_type):
+        if all(ct == "cartesian" for ct in coord_type):
             array = self.construct_array_cartesian(**kwargs)
-        elif all(type == "spherical" for type in coord_type):
+        elif all(ct == "spherical" for ct in coord_type):
             array = self.construct_array_spherical(**kwargs)
         elif isinstance(coord_type, (list, tuple)):
             array = self.construct_array_mix(coord_type, **kwargs)
diff --git a/gbasis/base_one.py b/gbasis/base_one.py
index e8270627..5da708c3 100644
--- a/gbasis/base_one.py
+++ b/gbasis/base_one.py
@@ -281,9 +281,9 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
             If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if all(type == "cartesian" for type in coord_type):
+        if all(ct == "cartesian" for ct in coord_type):
             array = self.construct_array_cartesian(**kwargs)
-        elif all(type == "spherical" for type in coord_type):
+        elif all(ct == "spherical" for ct in coord_type):
             array = self.construct_array_spherical(**kwargs)
         elif isinstance(coord_type, (list, tuple)):
             array = self.construct_array_mix(coord_type, **kwargs)
diff --git a/gbasis/base_two_asymm.py b/gbasis/base_two_asymm.py
index f1a86261..8cb18f59 100644
--- a/gbasis/base_two_asymm.py
+++ b/gbasis/base_two_asymm.py
@@ -421,9 +421,13 @@ def construct_array_lincomb(
             If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if all(type_one == "cartesian" for type_one in coord_type_one) and all(type_two == "cartesian" for type_two in coord_type_two):
+        if all(ct_one == "cartesian" for ct_one in coord_type_one) and all(
+            ct_two == "cartesian" for ct_two in coord_type_two
+        ):
             array = self.construct_array_cartesian(**kwargs)
-        elif all(type_one == "spherical" for type_one in coord_type_one) and all(type_two == "spherical" for type_two in coord_type_two):
+        elif all(ct_one == "spherical" for ct_one in coord_type_one) and all(
+            ct_two == "spherical" for ct_two in coord_type_two
+        ):
             array = self.construct_array_spherical(**kwargs)
         else:
             if coord_type_one in ["cartesian", "spherical"]:
diff --git a/gbasis/base_two_symm.py b/gbasis/base_two_symm.py
index 617914cb..a373d133 100644
--- a/gbasis/base_two_symm.py
+++ b/gbasis/base_two_symm.py
@@ -391,9 +391,9 @@ def construct_array_lincomb(self, transform, coord_type, **kwargs):
             If `coord_type` is not a list/tuple of the strings 'cartesian' or 'spherical'.
 
         """
-        if all(type == "cartesian" for type in coord_type):
+        if all(ct == "cartesian" for ct in coord_type):
             array = self.construct_array_cartesian(**kwargs)
-        elif all(type == "spherical" for type in coord_type):
+        elif all(ct == "spherical" for ct in coord_type):
             array = self.construct_array_spherical(**kwargs)
         elif isinstance(coord_type, (list, tuple)):
             array = self.construct_array_mix(coord_type, **kwargs)

From 610dd8f8d9777244ba7a7836293cf9762793a30e Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 10 Jan 2024 17:22:44 -0500
Subject: [PATCH 13/17] Allow 'c' or 'p' for make_contractions coord_types

---
 gbasis/parsers.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/gbasis/parsers.py b/gbasis/parsers.py
index 17727d28..5b0c5bb0 100644
--- a/gbasis/parsers.py
+++ b/gbasis/parsers.py
@@ -176,11 +176,11 @@ def make_contractions(basis_dict, atoms, coords, coord_types='spherical'):
         Atoms at which the contractions are centered.
     coords : np.ndarray(N, 3)
         Coordinates of each atom.
-    coord_type : {"cartesian", list/tuple of "cartesian" or "spherical", "spherical"}
+    coord_type : {"cartesian"/"c", list/tuple of "cartesian"/"c" or "spherical"/"p", "spherical"/"p"}
         Types of the coordinate system for the contractions.
-        If "cartesian", then all of the contractions are treated as Cartesian contractions.
-        If "spherical", then all of the contractions are treated as spherical contractions.
-        If list/tuple, then each entry must be a "cartesian" or "spherical" to specify the
+        If "cartesian" or "c", then all of the contractions are treated as Cartesian contractions.
+        If "spherical" or "p", then all of the contractions are treated as spherical contractions.
+        If list/tuple, then each entry must be a "cartesian" (or "c") or "spherical" (or "p") to specify the
         coordinate type of each `GeneralizedContractionShell` instance.
         Default value is "spherical".
 
@@ -214,12 +214,12 @@ def make_contractions(basis_dict, atoms, coords, coord_types='spherical'):
         for angmom, exps, coeffs in basis_dict[atom]:
             if type(coord_types) == str:
                 # if coord_types given as a single string, assign the specified type to all contractions for all atoms
-                if coord_types == "spherical":
+                if coord_types == "spherical" or coord_types == "p":
                     basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, 'p'))
-                elif coord_types == "cartesian":
+                elif coord_types == "cartesian" or coord_types == "c":
                     basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, 'c'))
                 else:
-                    raise ValueError("If coord_types is a string, it must be either 'spherical' or 'cartesian'.")
+                    raise ValueError("If coord_types is a string, it must be either 'spherical'/'p' or 'cartesian'/'c'.")
             elif type(coord_types) == list:
                 # if coord_types given as a list, assign the specified type to each atom's contractions individually
                 if len_coord_types == sum([len(basis_dict[i]) for i in atoms]):

From 4dedd3df64330c0f74be0580b65c4a210198913c Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 10 Jan 2024 18:36:41 -0500
Subject: [PATCH 14/17] Make coord_type arg required

---
 gbasis/contractions.py | 4 ++--
 gbasis/parsers.py      | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/gbasis/contractions.py b/gbasis/contractions.py
index 42a34564..1252987c 100644
--- a/gbasis/contractions.py
+++ b/gbasis/contractions.py
@@ -116,7 +116,7 @@ class GeneralizedContractionShell:
 
     """
 
-    def __init__(self, angmom, coord, coeffs, exps, coord_type="p"):
+    def __init__(self, angmom, coord, coeffs, exps, coord_type):
         r"""Initialize a GeneralizedContractionShell instance.
 
         Parameters
@@ -137,7 +137,7 @@ def __init__(self, angmom, coord, coeffs, exps, coord_type="p"):
             dimension.
         exps : np.ndarray(K,)
             Exponents of the primitives, :math:`\{\alpha_i\}_{i=1}^K`.
-        coord_type : str, optional
+        coord_type : str
             Coordinate type of the contraction. Options include "cartesian" or "c" and
             "spherical" or "p".
 
diff --git a/gbasis/parsers.py b/gbasis/parsers.py
index 5b0c5bb0..636670d3 100644
--- a/gbasis/parsers.py
+++ b/gbasis/parsers.py
@@ -165,7 +165,7 @@ def parse_gbs(gbs_basis_file):
     return output
 
 
-def make_contractions(basis_dict, atoms, coords, coord_types='spherical'):
+def make_contractions(basis_dict, atoms, coords, coord_types):
     """Return the contractions that correspond to the given atoms for the given basis.
 
     Parameters

From 1e8931c13218d9f78bb1d202ac17dc362ccffd7f Mon Sep 17 00:00:00 2001
From: maximilianvz <max.p.vanzyl@gmail.com>
Date: Wed, 10 Jan 2024 18:45:26 -0500
Subject: [PATCH 15/17] Adjust tests for required coord_type arg

---
 tests/test_angular_momentum.py            |  4 +--
 tests/test_base.py                        | 10 +++---
 tests/test_base_four_symm.py              | 34 +++++++++----------
 tests/test_base_one.py                    | 20 ++++++------
 tests/test_base_two_asymm.py              | 30 ++++++++---------
 tests/test_base_two_symm.py               | 40 +++++++++++------------
 tests/test_contractions.py                | 11 ++++---
 tests/test_density.py                     | 32 +++++++++---------
 tests/test_electron_repulsion.py          | 14 ++++----
 tests/test_eval.py                        |  4 +--
 tests/test_eval_deriv.py                  |  6 ++--
 tests/test_kinetic_energy.py              |  4 +--
 tests/test_moment.py                      |  6 ++--
 tests/test_momentum.py                    |  4 +--
 tests/test_nuclear_electron_attraction.py |  2 +-
 tests/test_one_elec_int.py                | 16 ++++-----
 tests/test_overlap.py                     | 14 ++++----
 tests/test_parsers.py                     | 14 ++++----
 tests/test_point_charge.py                | 12 +++----
 tests/test_stress_tensor.py               |  8 ++---
 tests/test_wrappers.py                    |  2 +-
 21 files changed, 144 insertions(+), 143 deletions(-)

diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py
index 9f73d81f..87d44a99 100644
--- a/tests/test_angular_momentum.py
+++ b/tests/test_angular_momentum.py
@@ -14,10 +14,10 @@
 def test_angular_momentum_construct_array_contraction():
     """Test integrals.angular_momentum.angular_momentumIntegral.construct_array_contraction."""
     test_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02])
+        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02]), 'spherical'
     )
 
     # copied the code it is testing
diff --git a/tests/test_base.py b/tests/test_base.py
index cf3c00a6..6a4e8d9f 100644
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -10,7 +10,7 @@ def test_init():
     """Test base.BaseGaussianRelatedArray."""
     Test = disable_abstract(BaseGaussianRelatedArray)  # noqa: N806
     test = skip_init(Test)
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     assert not hasattr(test, "_axes_contractions")
     with pytest.raises(TypeError):
         Test.__init__(test, set([contractions]))
@@ -38,7 +38,7 @@ def test_contruct_array_contraction():
             "construct_array_contraction": BaseGaussianRelatedArray.construct_array_contraction
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
@@ -52,7 +52,7 @@ def test_contruct_array_cartesian():
             "construct_array_cartesian": BaseGaussianRelatedArray.construct_array_cartesian
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
@@ -66,7 +66,7 @@ def test_contruct_array_spherical():
             "construct_array_spherical": BaseGaussianRelatedArray.construct_array_spherical
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
@@ -80,6 +80,6 @@ def test_contruct_array_lincomb():
             "construct_array_lincomb": BaseGaussianRelatedArray.construct_array_lincomb
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
diff --git a/tests/test_base_four_symm.py b/tests/test_base_four_symm.py
index 381e18bd..92d2548c 100644
--- a/tests/test_base_four_symm.py
+++ b/tests/test_base_four_symm.py
@@ -11,7 +11,7 @@ def test_init():
     """Test BaseFourIndexSymmetric.__init__."""
     Test = disable_abstract(BaseFourIndexSymmetric)  # noqa: N806
     test = skip_init(Test)
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test.__init__(test, [contractions])
     assert test._axes_contractions[0][0] == contractions
     with pytest.raises(TypeError):
@@ -21,7 +21,7 @@ def test_init():
 def test_contractions():
     """Test BaseFourIndexSymmetric.contractions."""
     Test = disable_abstract(BaseFourIndexSymmetric)  # noqa: N806
-    cont = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     test = Test([cont])
     assert test.contractions[0] == cont
 
@@ -35,15 +35,15 @@ def test_construct_array_contraction():
             "construct_array_contraction": BaseFourIndexSymmetric.construct_array_contraction
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
 
 def test_construct_array_cartesian():
     """Test BaseFourIndexSymmetric.construct_array_cartesian."""
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones((1, 1)), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([2, 3, 4]), np.ones((1, 1)), 2 * np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones((1, 1)), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([2, 3, 4]), np.ones((1, 1)), 2 * np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseFourIndexSymmetric,
         dict_overwrite={
@@ -91,7 +91,7 @@ def test_construct_array_cartesian():
 
 def test_construct_array_spherical():
     """Test BaseFourIndexSymmetric.construct_array_spherical."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -127,8 +127,8 @@ def test_construct_array_spherical():
     with pytest.raises(TypeError):
         test.construct_array_spherical(bad_keyword=3)
 
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     transform_one = generate_transformation(
         1, cont_one.angmom_components_cart, cont_one.angmom_components_sph, "left"
     )
@@ -400,7 +400,7 @@ def test_construct_array_spherical():
 
 def test_construct_array_mix():
     """Test BaseFourIndex.construct_array_mix."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     Test = disable_abstract(  # noqa: N806
         BaseFourIndexSymmetric,
@@ -423,8 +423,8 @@ def test_construct_array_mix():
         test.construct_array_cartesian(a=3), test.construct_array_mix(["cartesian"], a=3)
     )
 
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     Test = disable_abstract(  # noqa: N806
         BaseFourIndexSymmetric,
@@ -566,8 +566,8 @@ def construct_array_cont(self, cont_one, cont_two, cont_three, cont_four):
         BaseFourIndexSymmetric,
         dict_overwrite={"construct_array_contraction": construct_array_cont},
     )
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     # Remove the dependence on norm constants.
     cont_one.norm_cont = np.ones((1, cont_one.num_cart))
@@ -656,7 +656,7 @@ def construct_array_cont(self, cont_one, cont_two, cont_three, cont_four):
 
 def test_construct_array_lincomb():
     """Test BaseFourIndexSymmetric.construct_array_lincomb."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     sph_transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -732,8 +732,8 @@ def test_construct_array_lincomb():
             )
         },
     )
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     cont_one.norm_cont = np.ones((1, cont_one.num_cart))
     cont_two.norm_cont = np.ones((1, cont_two.num_cart))
     test = Test([cont_one, cont_two])
@@ -1250,7 +1250,7 @@ def angmom_components_sph(self):
             """Raise error in case undefined conventions are accessed."""
             raise NotImplementedError
 
-    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1))
+    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseFourIndexSymmetric,
         dict_overwrite={
diff --git a/tests/test_base_one.py b/tests/test_base_one.py
index 043166b7..df6b669d 100644
--- a/tests/test_base_one.py
+++ b/tests/test_base_one.py
@@ -11,7 +11,7 @@ def test_init():
     """Test BaseOneIndex.__init__."""
     Test = disable_abstract(BaseOneIndex)  # noqa: N806
     test = skip_init(Test)
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test.__init__(test, [contractions])
     assert test._axes_contractions == ((contractions,),)
     with pytest.raises(TypeError):
@@ -21,7 +21,7 @@ def test_init():
 def test_contractions():
     """Test BaseOneIndex.constractions."""
     Test = disable_abstract(BaseOneIndex)  # noqa: N806
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     test = Test([contractions])
     assert test.contractions[0] == contractions
 
@@ -33,14 +33,14 @@ def test_contruct_array_contraction():
         BaseOneIndex,
         dict_overwrite={"construct_array_contraction": BaseOneIndex.construct_array_contraction},
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
 
 def test_contruct_array_cartesian():
     """Test BaseOneIndex.construct_array_cartesian."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     contractions.norm_cont = np.ones((1, 5))
     Test = disable_abstract(  # noqa: N806
         BaseOneIndex,
@@ -78,7 +78,7 @@ def test_contruct_array_cartesian():
 
 def test_contruct_array_spherical():
     """Test BaseOneIndex.construct_array_spherical."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -108,7 +108,7 @@ def test_contruct_array_spherical():
         np.vstack([transform.dot(np.arange(9).reshape(3, 3)) * 2] * 2),
     )
 
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseOneIndex,
         dict_overwrite={
@@ -155,7 +155,7 @@ def test_contruct_array_spherical():
 
 def test_contruct_array_mix():
     """Test BaseOneIndex.construct_array_mix."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     Test = disable_abstract(  # noqa: N806
         BaseOneIndex,
@@ -190,7 +190,7 @@ def test_contruct_array_mix():
         test.construct_array_cartesian(a=3), test.construct_array_mix(["cartesian"] * 2, a=3)
     )
 
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseOneIndex,
         dict_overwrite={
@@ -243,7 +243,7 @@ def angmom_components_sph(self):
             """Raise error in case undefined conventions are accessed."""
             raise NotImplementedError
 
-    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1))
+    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseOneIndex,
         dict_overwrite={
@@ -260,7 +260,7 @@ def angmom_components_sph(self):
 
 def test_contruct_array_lincomb():
     """Test BaseOneIndex.construct_array_lincomb."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     sph_transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
diff --git a/tests/test_base_two_asymm.py b/tests/test_base_two_asymm.py
index c9dba68d..76910628 100644
--- a/tests/test_base_two_asymm.py
+++ b/tests/test_base_two_asymm.py
@@ -11,7 +11,7 @@ def test_init():
     """Test BaseTwoIndexAsymmetric.__init__."""
     Test = disable_abstract(BaseTwoIndexAsymmetric)  # noqa: N806
     test = skip_init(Test)
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test.__init__(test, [contractions], [contractions])
     assert test._axes_contractions == ((contractions,), (contractions,))
     with pytest.raises(TypeError):
@@ -23,8 +23,8 @@ def test_init():
 def test_contractions_one():
     """Test BaseTwoIndexAsymmetric.constractions_one."""
     Test = disable_abstract(BaseTwoIndexAsymmetric)  # noqa: N806
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     test = Test([cont_one], [cont_two])
     assert test.contractions_one[0] == cont_one
 
@@ -32,8 +32,8 @@ def test_contractions_one():
 def test_contractions_two():
     """Test BaseTwoIndexAsymmetric.constractions_two."""
     Test = disable_abstract(BaseTwoIndexAsymmetric)  # noqa: N806
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     test = Test([cont_one], [cont_two])
     assert test.contractions_two[0] == cont_two
 
@@ -47,14 +47,14 @@ def test_contruct_array_contraction():
             "construct_array_contraction": BaseTwoIndexAsymmetric.construct_array_contraction
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
 
 def test_contruct_array_cartesian():
     """Test BaseTwoIndexAsymmetric.construct_array_cartesian."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexAsymmetric,
         dict_overwrite={
@@ -80,8 +80,8 @@ def test_contruct_array_cartesian():
             )
         },
     )
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     cont_one.norm_cont = np.ones((1, 2))
     cont_two.norm_cont = np.ones((1, 5))
     test = Test([cont_one, cont_one], [cont_two])
@@ -105,7 +105,7 @@ def test_contruct_array_cartesian():
 
 def test_contruct_array_spherical():
     """Test BaseTwoIndexAsymmetric.construct_array_spherical."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -194,7 +194,7 @@ def test_contruct_array_spherical():
 
 def test_contruct_array_mix():
     """Test BaseTwoIndexAsymmetric.construct_array_mix."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexAsymmetric,
@@ -347,8 +347,8 @@ def construct_array_cont(self, cont_one, cont_two):
         BaseTwoIndexAsymmetric,
         dict_overwrite={"construct_array_contraction": construct_array_cont},
     )
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     # Remove the dependence on norm constants.
     cont_one.norm_cont = np.ones((1, cont_one.num_cart))
@@ -393,7 +393,7 @@ def construct_array_cont(self, cont_one, cont_two):
 
 def test_contruct_array_lincomb():
     """Test BaseTwoIndexAsymmetric.construct_array_lincomb."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     sph_transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -556,7 +556,7 @@ def angmom_components_sph(self):
             """Raise error in case undefined conventions are accessed."""
             raise NotImplementedError
 
-    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1))
+    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexAsymmetric,
         dict_overwrite={
diff --git a/tests/test_base_two_symm.py b/tests/test_base_two_symm.py
index a5add74d..5fcecfaf 100644
--- a/tests/test_base_two_symm.py
+++ b/tests/test_base_two_symm.py
@@ -12,7 +12,7 @@ def test_init():
     """Test BaseTwoIndexSymmetric.__init__."""
     Test = disable_abstract(BaseTwoIndexSymmetric)  # noqa: N806
     test = skip_init(Test)
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test.__init__(test, [contractions])
     assert test._axes_contractions[0][0] == contractions
     with pytest.raises(TypeError):
@@ -22,7 +22,7 @@ def test_init():
 def test_contractions():
     """Test BaseTwoIndexSymmetric.constractions."""
     Test = disable_abstract(BaseTwoIndexSymmetric)  # noqa: N806
-    cont = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     test = Test([cont])
     assert test.contractions[0] == cont
 
@@ -36,7 +36,7 @@ def test_contruct_array_contraction():
             "construct_array_contraction": BaseTwoIndexSymmetric.construct_array_contraction
         },
     )
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     with pytest.raises(TypeError):
         Test([contractions])
 
@@ -49,7 +49,7 @@ def test_contruct_array_contraction():
 # to the tril blocks because the ordering is different.
 def test_contruct_array_cartesian():
     """Test BaseTwoIndexSymmetric.construct_array_cartesian."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexSymmetric,
         dict_overwrite={
@@ -67,8 +67,8 @@ def test_contruct_array_cartesian():
     assert np.allclose(test.construct_array_cartesian(), np.ones((4, 4)) * 2)
     assert np.allclose(test.construct_array_cartesian(a=3), np.ones((4, 4)) * 3)
 
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexSymmetric,
         dict_overwrite={
@@ -229,7 +229,7 @@ def test_contruct_array_cartesian():
 # to the tril blocks because the ordering is different.
 def test_contruct_array_spherical():
     """Test BaseTwoIndexSymmetric.construct_array_spherical."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -256,8 +256,8 @@ def test_contruct_array_spherical():
     with pytest.raises(TypeError):
         test.construct_array_spherical(bad_keyword=3)
 
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     transform_one = generate_transformation(
         1, cont_one.angmom_components_cart, cont_one.angmom_components_sph, "left"
     )
@@ -400,7 +400,7 @@ def test_contruct_array_spherical():
 
 def test_contruct_array_mix():
     """Test BaseTwoIndex.construct_array_mix."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexSymmetric,
@@ -421,8 +421,8 @@ def test_contruct_array_mix():
         test.construct_array_cartesian(a=3), test.construct_array_mix(["cartesian"], a=3)
     )
 
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexSymmetric,
@@ -489,7 +489,7 @@ def test_contruct_array_mix():
 
 def test_contruct_array_lincomb():
     """Test BaseTwoIndexSymmetric.construct_array_lincomb."""
-    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    contractions = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     sph_transform = generate_transformation(
         1, contractions.angmom_components_cart, contractions.angmom_components_sph, "left"
     )
@@ -548,8 +548,8 @@ def test_contruct_array_lincomb():
             )
         },
     )
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     cont_one.norm_cont = np.ones((1, cont_one.num_cart))
     cont_two.norm_cont = np.ones((1, cont_two.num_cart))
     test = Test([cont_one, cont_two])
@@ -733,8 +733,8 @@ def construct_array_cont(self, cont_one, cont_two):
         BaseTwoIndexSymmetric,
         dict_overwrite={"construct_array_contraction": construct_array_cont},
     )
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
 
     # Remove the dependence on norm constants.
     cont_one.norm_cont = np.ones((1, cont_one.num_cart))
@@ -777,8 +777,8 @@ def construct_array_cont(self, cont_one, cont_two):
 
 def test_compare_two_asymm():
     """Test BaseTwoIndexSymmetric by comparing it against BaseTwoIndexAsymmetric."""
-    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1))
-    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1))
+    cont_one = GeneralizedContractionShell(1, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
+    cont_two = GeneralizedContractionShell(2, np.array([1, 2, 3]), np.ones(1), np.ones(1), 'spherical')
     sph_orb_transform = np.random.rand(8, 8)
     cart_orb_transform = np.random.rand(9, 9)
 
@@ -852,7 +852,7 @@ def angmom_components_sph(self):
             """Raise error in case undefined conventions are accessed."""
             raise NotImplementedError
 
-    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1))
+    contractions = SpecialShell(1, np.array([1, 2, 3]), np.ones((1, 2)), np.ones(1), 'spherical')
     Test = disable_abstract(  # noqa: N806
         BaseTwoIndexSymmetric,
         dict_overwrite={
diff --git a/tests/test_contractions.py b/tests/test_contractions.py
index 924e3d57..198e2d49 100644
--- a/tests/test_contractions.py
+++ b/tests/test_contractions.py
@@ -169,6 +169,7 @@ def tests_init():
         np.array([0, 1, 2]),
         np.array([1, 2, 3, 4], dtype=float),
         np.array([5, 6, 7, 8], dtype=float),
+        'spherical'
     )
     assert test._angmom == 1
     assert np.allclose(test._coord, np.array([0, 1, 2]))
@@ -225,9 +226,9 @@ def test_angmom_components_sph():
 # TODO: add more tests
 def test_norm_prim_cart():
     """Test GeneralizedContractionShell.norm_prim_cart."""
-    test = GeneralizedContractionShell(0, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]))
+    test = GeneralizedContractionShell(0, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]), 'spherical')
     assert np.isclose(test.norm_prim_cart, 0.2519794355383807303479140)
-    test = GeneralizedContractionShell(3, np.array([0, 0, 0]), np.array([1.0]), np.array([0.5]))
+    test = GeneralizedContractionShell(3, np.array([0, 0, 0]), np.array([1.0]), np.array([0.5]), 'spherical')
     assert np.isclose(test.norm_prim_cart[7], 0.6920252830162908851679097)
 
 
@@ -260,14 +261,14 @@ def test_num_seg_cont():
 
 def test_assign_norm_cont():
     """Test GeneralizedContractionShell.assign_norm_cont."""
-    test = GeneralizedContractionShell(0, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]))
+    test = GeneralizedContractionShell(0, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]), 'spherical')
     test.assign_norm_cont()
     assert np.allclose(test.norm_cont, 1)
 
-    test = GeneralizedContractionShell(1, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]))
+    test = GeneralizedContractionShell(1, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]), 'spherical')
     test.assign_norm_cont()
     assert np.allclose(test.norm_cont, 1)
 
-    test = GeneralizedContractionShell(2, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]))
+    test = GeneralizedContractionShell(2, np.array([0, 0, 0]), np.array([1.0]), np.array([0.25]), 'spherical')
     test.assign_norm_cont()
     assert np.allclose(test.norm_cont, 1)
diff --git a/tests/test_density.py b/tests/test_density.py
index 9f61c51f..301d6225 100644
--- a/tests/test_density.py
+++ b/tests/test_density.py
@@ -63,7 +63,7 @@ def test_evaluate_density_using_evaluated_orbs():
 def test_evaluate_density():
     """Test gbasis.evals.density.evaluate_density."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     transform = np.random.rand(14, 18)
     density = np.random.rand(14, 14)
     density += density.T
@@ -79,7 +79,7 @@ def test_evaluate_density():
 def test_evaluate_deriv_density():
     """Test gbasis.evals.density.evaluate_deriv_density."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     transform = np.random.rand(14, 18)
     density = np.random.rand(14, 14)
     density += density.T
@@ -223,7 +223,7 @@ def test_evaluate_deriv_density():
 def test_evaluate_density_gradient():
     """Test gbasis.evals.density.evaluate_density_gradient."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     transform = np.random.rand(14, 18)
     density = np.random.rand(14, 14)
     density += density.T
@@ -283,8 +283,8 @@ def test_evaluate_density_horton():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     points = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_density = np.load(find_datafile("data_horton_hhe_sph_density.npy"))
 
@@ -306,8 +306,8 @@ def test_evaluate_density_gradient_horton():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     points = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_density_gradient = np.load(find_datafile("data_horton_hhe_sph_density_gradient.npy"))
 
@@ -330,8 +330,8 @@ def test_evaluate_hessian_deriv_horton():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     points = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_density_hessian = np.zeros((10 ** 3, 3, 3))
     horton_density_hessian[:, [0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]] = np.load(
@@ -360,8 +360,8 @@ def test_evaluate_laplacian_deriv_horton():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     points = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_density_laplacian = np.load(find_datafile("data_horton_hhe_sph_density_laplacian.npy"))
 
@@ -384,8 +384,8 @@ def test_evaluate_posdef_kinetic_energy_density():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     points = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_density_kinetic_density = np.load(
         find_datafile("data_horton_hhe_sph_posdef_kinetic_density.npy")
@@ -412,8 +412,8 @@ def test_evaluate_general_kinetic_energy_density_horton():
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     # NOTE: used HORTON's conversion factor for angstroms to bohr
     points = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     horton_density_kinetic_density = np.load(
         find_datafile("data_horton_hhe_sph_posdef_kinetic_density.npy")
@@ -435,7 +435,7 @@ def test_evaluate_general_kinetic_energy_density():
     """Test density.evaluate_general_kinetic_energy_density."""
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     points = np.array([[0, 0, 0]])
-    basis = make_contractions(basis_dict, ["H"], points)
+    basis = make_contractions(basis_dict, ["H"], points, 'spherical')
     points = np.random.rand(10, 3)
 
     with pytest.raises(TypeError):
diff --git a/tests/test_electron_repulsion.py b/tests/test_electron_repulsion.py
index c76c73b0..44ab0dd4 100644
--- a/tests/test_electron_repulsion.py
+++ b/tests/test_electron_repulsion.py
@@ -17,13 +17,13 @@
 def test_construct_array_contraction():
     """Test integrals.electron_repulsion.ElectronRepulsionIntegral.construct_array_contraction."""
     coord_one = np.array([0.5, 1, 1.5])
-    cont_one = GeneralizedContractionShell(0, coord_one, np.array([1.0]), np.array([0.1]))
+    cont_one = GeneralizedContractionShell(0, coord_one, np.array([1.0]), np.array([0.1]), 'spherical')
     coord_two = np.array([1.5, 2, 3])
-    cont_two = GeneralizedContractionShell(0, coord_two, np.array([3.0]), np.array([0.2]))
+    cont_two = GeneralizedContractionShell(0, coord_two, np.array([3.0]), np.array([0.2]), 'spherical')
     coord_three = np.array([2.5, 3, 4])
-    cont_three = GeneralizedContractionShell(0, coord_three, np.array([3.0]), np.array([0.2]))
+    cont_three = GeneralizedContractionShell(0, coord_three, np.array([3.0]), np.array([0.2]), 'spherical')
     coord_four = np.array([3.5, 4, 5])
-    cont_four = GeneralizedContractionShell(0, coord_four, np.array([3.0]), np.array([0.2]))
+    cont_four = GeneralizedContractionShell(0, coord_four, np.array([3.0]), np.array([0.2]), 'spherical')
 
     with pytest.raises(TypeError):
         ElectronRepulsionIntegral.construct_array_contraction(None, cont_two, cont_three, cont_four)
@@ -123,10 +123,10 @@ def test_electron_repulsion_cartesian_horton_custom_hhe():
     basis = make_contractions(basis_dict, ["H", "He"], coords, 'cartesian')
     basis = [HortonContractions(i.angmom, i.coord, i.coeffs[:, 0], i.exps, i.coord_type) for i in basis[:8]]
     basis[0] = HortonContractions(
-        basis[0].angmom, basis[0].coord, basis[0].coeffs[3:], basis[0].exps[3:]
+        basis[0].angmom, basis[0].coord, basis[0].coeffs[3:], basis[0].exps[3:], basis[0].coord_type
     )
     basis[4] = HortonContractions(
-        basis[4].angmom, basis[4].coord, basis[4].coeffs[4:], basis[4].exps[4:]
+        basis[4].angmom, basis[4].coord, basis[4].coeffs[4:], basis[4].exps[4:], basis[4].coord_type
     )
     basis.pop(3)
     basis.pop(2)
@@ -194,7 +194,7 @@ def test_electron_repulsion_mix():
 def test_electron_repulsion_lincomb():
     """Test gbasis.integrals.electron_repulsion.electron_repulsion_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]), "spherical")
 
     erep_obj = ElectronRepulsionIntegral(basis)
     transform = np.random.rand(3, 5)
diff --git a/tests/test_eval.py b/tests/test_eval.py
index 954477d1..640a6d6b 100644
--- a/tests/test_eval.py
+++ b/tests/test_eval.py
@@ -12,7 +12,7 @@
 def test_evaluate_construct_array_contraction():
     """Test gbasis.evals.eval.Eval.construct_array_contraction."""
     test = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     answer = np.array(
         [
@@ -141,7 +141,7 @@ def test_evaluate_basis_horton():
     """Test gbasis.evals.eval.evaluate_basis against horton results."""
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     points = np.array([[0, 0, 0], [0.8, 0, 0]])
-    basis = make_contractions(basis_dict, ["H", "He"], points)
+    basis = make_contractions(basis_dict, ["H", "He"], points, 'spherical')
     cartesian_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, 'cartesian') for i in basis]
     spherical_basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, 'spherical') for i in basis]
 
diff --git a/tests/test_eval_deriv.py b/tests/test_eval_deriv.py
index e00c07c1..c32f6f9c 100644
--- a/tests/test_eval_deriv.py
+++ b/tests/test_eval_deriv.py
@@ -16,7 +16,7 @@ def test_evaluate_deriv_construct_array_contraction():
     points = np.array([[2, 3, 4]])
     orders = np.array([0, 0, 0])
     contractions = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     with pytest.raises(TypeError):
         EvalDeriv.construct_array_contraction(
@@ -53,7 +53,7 @@ def test_evaluate_deriv_construct_array_contraction():
         orders[k] = 1
 
         test = GeneralizedContractionShell(
-            1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+            1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
         )
         answer = np.array(
             [
@@ -93,7 +93,7 @@ def test_evaluate_deriv_construct_array_contraction():
         orders[l] += 1
 
         test = GeneralizedContractionShell(
-            1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+            1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
         )
         answer = np.array(
             [
diff --git a/tests/test_kinetic_energy.py b/tests/test_kinetic_energy.py
index 7bad067c..286dc2ec 100644
--- a/tests/test_kinetic_energy.py
+++ b/tests/test_kinetic_energy.py
@@ -12,10 +12,10 @@
 def test_kinetic_energy_construct_array_contraction():
     """Test gbasis.integrals.kinetic_energy.KineticEnergyIntegral.construct_array_contraction."""
     test_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02])
+        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02]), 'spherical'
     )
     answer = np.array(
         [
diff --git a/tests/test_moment.py b/tests/test_moment.py
index 2195f828..b191e14e 100644
--- a/tests/test_moment.py
+++ b/tests/test_moment.py
@@ -12,10 +12,10 @@
 def test_moment_construct_array_contraction():
     """Test gbasis.integrals.moment.Moment.construct_array_contraction."""
     test_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02])
+        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02]), 'spherical'
     )
     answer = np.array(
         [
@@ -270,7 +270,7 @@ def test_moment_mix():
 def test_moment_spherical_lincomb():
     """Test gbasis.integrals.moment.moment_spherical_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     moment_obj = Moment(basis)
     transform = np.random.rand(14, 18)
     assert np.allclose(
diff --git a/tests/test_momentum.py b/tests/test_momentum.py
index 5bd2c865..645e6c63 100644
--- a/tests/test_momentum.py
+++ b/tests/test_momentum.py
@@ -12,10 +12,10 @@
 def test_momentum_construct_array_contraction():
     """Test gbasis.integrals.momentum.MomentumIntegral.construct_array_contraction."""
     test_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02])
+        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02]), 'spherical'
     )
     test = MomentumIntegral.construct_array_contraction(test_one, test_two).squeeze()
     answer = np.array(
diff --git a/tests/test_nuclear_electron_attraction.py b/tests/test_nuclear_electron_attraction.py
index 2bbf00aa..61c29de0 100644
--- a/tests/test_nuclear_electron_attraction.py
+++ b/tests/test_nuclear_electron_attraction.py
@@ -101,7 +101,7 @@ def test_nuclear_electron_attraction_mix():
 def test_nuclear_electron_attraction_lincomb():
     """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_lincomb."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]),)
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
 
     nuclear_coords = np.random.rand(5, 3)
     nuclear_charges = np.random.rand(5)
diff --git a/tests/test_one_elec_int.py b/tests/test_one_elec_int.py
index f9a15a05..244f52a4 100644
--- a/tests/test_one_elec_int.py
+++ b/tests/test_one_elec_int.py
@@ -46,10 +46,10 @@ def boys_func(order, weighted_dist):
 def test_compute_one_elec_int_v_recursion():
     """Test vertical recursion in _one_elec_int._compute_one_elec_integrals."""
     contr_one = GeneralizedContractionShell(
-        3, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.25])
+        3, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.25]), 'spherical'
     )
     contr_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.02, 0.01])
+        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.02, 0.01]), 'spherical'
     )
     coord_a = contr_one.coord
     angmom_a = contr_one.angmom
@@ -98,10 +98,10 @@ def test_compute_one_elec_int_s_type():
     """Test _one_elec_int._compute_one_electron_integrals for s-type primitives."""
     # GeneralizedContractionShell(angmom, coord, charge, coeffs, exps)
     s_type_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0]), np.array([0.1])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0]), np.array([0.1]), 'spherical'
     )
     s_type_two = GeneralizedContractionShell(
-        1, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.02])
+        1, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.02]), 'spherical'
     )
     coord_a = s_type_one.coord
     angmom_a = s_type_one.angmom
@@ -182,10 +182,10 @@ def test_compute_one_elec_int_multiple_contractions():
     """Test _one_elec_int._compute_one_electron_integrals for s-type contractions."""
     # GeneralizedContractionShell(angmom, coord, charge, coeffs, exps)
     contr_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.25])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.25]), 'spherical'
     )
     contr_two = GeneralizedContractionShell(
-        1, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.02, 0.01])
+        1, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.02, 0.01]), 'spherical'
     )
     coord_a = contr_one.coord
     angmom_a = contr_one.angmom
@@ -227,10 +227,10 @@ def test_compute_one_elec_int_multiple_contractions():
 def test_compute_one_elec_int_generalized_contraction():
     """Test _one_elec_int._compute_one_electron_integrals for generalized contractions."""
     contr_one = GeneralizedContractionShell(
-        3, np.array([0.5, 1, 1.5]), np.array([[1.0, 2.0], [1.5, 2.5]]), np.array([0.1, 0.25])
+        3, np.array([0.5, 1, 1.5]), np.array([[1.0, 2.0], [1.5, 2.5]]), np.array([0.1, 0.25]), 'spherical'
     )
     contr_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([[3.0, 4.0], [3.5, 4.5]]), np.array([0.02, 0.01])
+        2, np.array([1.5, 2, 3]), np.array([[3.0, 4.0], [3.5, 4.5]]), np.array([0.02, 0.01]), 'spherical'
     )
     coord_a = contr_one.coord
     angmom_a = contr_one.angmom
diff --git a/tests/test_overlap.py b/tests/test_overlap.py
index 8a51fcf9..ac84b232 100644
--- a/tests/test_overlap.py
+++ b/tests/test_overlap.py
@@ -13,10 +13,10 @@
 def test_overlap_construct_array_contraction():
     """Test gbasis.integrals.overlap.Overlap.construct_array_contraction."""
     test_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02])
+        2, np.array([1.5, 2, 3]), np.array([3.0, 4.0]), np.array([0.2, 0.02]), 'spherical'
     )
     answer = np.array(
         [
@@ -125,7 +125,7 @@ def test_overlap_cartesian_norm_anorcc():
     """
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     overlap_obj = Overlap(basis)
     assert np.allclose(np.diag(overlap_obj.construct_array_cartesian()), 1)
 
@@ -139,7 +139,7 @@ def test_overlap_spherical_norm_sto6g():
     """
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
     overlap_obj = Overlap(basis)
     assert np.allclose(np.diag(overlap_obj.construct_array_spherical()), 1)
 
@@ -152,11 +152,11 @@ def test_overlap_spherical_norm_anorcc():
     """
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
 
-    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["C"], np.array([[0, 0, 0]]), 'cartesian')
     overlap_obj = Overlap(basis)
     assert np.allclose(np.diag(overlap_obj.construct_array_cartesian()), 1)
 
-    basis = make_contractions(basis_dict, ["Xe"], np.array([[0, 0, 0]]))
+    basis = make_contractions(basis_dict, ["Xe"], np.array([[0, 0, 0]]), 'cartesian')
     overlap_obj = Overlap(basis)
     assert np.allclose(np.diag(overlap_obj.construct_array_cartesian()), 1)
 
@@ -170,7 +170,7 @@ def test_overlap_cartesian_norm_sto6g():
     """
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
 
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'spherical')
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), 'cartesian')
     overlap_obj = Overlap(basis)
     assert np.allclose(np.diag(overlap_obj.construct_array_cartesian()), 1)
 
diff --git a/tests/test_parsers.py b/tests/test_parsers.py
index b2cde802..10059787 100644
--- a/tests/test_parsers.py
+++ b/tests/test_parsers.py
@@ -762,24 +762,24 @@ def test_make_contractions():
     """Test gbasis.contractions.make_contractions."""
     basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
     with pytest.raises(TypeError):
-        make_contractions(basis_dict, {"H", "H"}, np.array([[0, 0, 0], [1, 1, 1]]))
+        make_contractions(basis_dict, {"H", "H"}, np.array([[0, 0, 0], [1, 1, 1]]), 'spherical')
     with pytest.raises(TypeError):
-        make_contractions(basis_dict, [0, 0], np.array([[0, 0, 0], [1, 1, 1]]))
+        make_contractions(basis_dict, [0, 0], np.array([[0, 0, 0], [1, 1, 1]]), 'spherical')
 
     with pytest.raises(TypeError):
-        make_contractions(basis_dict, ["H", "H"], [[0, 0, 0], [1, 1, 1]])
+        make_contractions(basis_dict, ["H", "H"], [[0, 0, 0], [1, 1, 1]], 'spherical')
     with pytest.raises(TypeError):
-        make_contractions(basis_dict, ["H", "H"], np.array([0, 0, 0, 1, 1, 1]))
+        make_contractions(basis_dict, ["H", "H"], np.array([0, 0, 0, 1, 1, 1]), 'spherical')
     with pytest.raises(TypeError):
-        make_contractions(basis_dict, ["H", "H"], np.array([[0, 0, 0, 2], [1, 1, 1, 2]]))
+        make_contractions(basis_dict, ["H", "H"], np.array([[0, 0, 0, 2], [1, 1, 1, 2]]), 'spherical')
 
     with pytest.raises(ValueError):
-        make_contractions(basis_dict, ["H", "H", "H"], np.array([[0, 0, 0], [1, 1, 1]]))
+        make_contractions(basis_dict, ["H", "H", "H"], np.array([[0, 0, 0], [1, 1, 1]]), 'spherical')
 
     with pytest.raises(TypeError):
         make_contractions(basis_dict, ["H", "H"], np.array([[0, 0, 0], [1, 1, 1]]), [0, 0])
 
-    test = make_contractions(basis_dict, ["H", "H"], np.array([[0, 0, 0], [1, 1, 1]]))
+    test = make_contractions(basis_dict, ["H", "H"], np.array([[0, 0, 0], [1, 1, 1]]), 'spherical')
     assert isinstance(test, tuple)
     assert len(test) == 2
     assert test[0].angmom == 0
diff --git a/tests/test_point_charge.py b/tests/test_point_charge.py
index 3add66ce..074656a9 100644
--- a/tests/test_point_charge.py
+++ b/tests/test_point_charge.py
@@ -33,9 +33,9 @@ def test_boys_func():
 def test_construct_array_contraction():
     """Test gbasis.integrals.point_charge.PointChargeIntegral.construct_array_contraction."""
     coord_one = np.array([0.5, 1, 1.5])
-    test_one = GeneralizedContractionShell(0, coord_one, np.array([1.0]), np.array([0.1]))
+    test_one = GeneralizedContractionShell(0, coord_one, np.array([1.0]), np.array([0.1]), 'spherical')
     coord_two = np.array([1.5, 2, 3])
-    test_two = GeneralizedContractionShell(0, coord_two, np.array([3.0]), np.array([0.2]))
+    test_two = GeneralizedContractionShell(0, coord_two, np.array([3.0]), np.array([0.2]), 'spherical')
     coord = np.array([[0, 0, 0]])
     charge = np.array([1])
     coord_wac = (0.1 * coord_one + 0.2 * coord_two) / (0.1 + 0.2)
@@ -70,10 +70,10 @@ def test_construct_array_contraction():
     )
 
     test_one = GeneralizedContractionShell(
-        1, np.array([0.5, 1, 1.5]), np.array([1.0]), np.array([0.1])
+        1, np.array([0.5, 1, 1.5]), np.array([1.0]), np.array([0.1]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        0, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.2])
+        0, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.2]), 'spherical'
     )
     v_000_000 = [
         2
@@ -97,10 +97,10 @@ def test_construct_array_contraction():
     )
 
     test_one = GeneralizedContractionShell(
-        0, np.array([0.5, 1, 1.5]), np.array([1.0]), np.array([0.1])
+        0, np.array([0.5, 1, 1.5]), np.array([1.0]), np.array([0.1]), 'spherical'
     )
     test_two = GeneralizedContractionShell(
-        1, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.2])
+        1, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.2]), 'spherical'
     )
     assert np.allclose(
         PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge).ravel(),
diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py
index 9a90ba0b..5aef2759 100644
--- a/tests/test_stress_tensor.py
+++ b/tests/test_stress_tensor.py
@@ -19,7 +19,7 @@ def test_evaluate_stress_tensor():
     """Test gbasis.evals.stress_tensor.evaluate_stress_tensor."""
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     coords = np.array([[0, 0, 0]])
-    basis = make_contractions(basis_dict, ["H"], coords)
+    basis = make_contractions(basis_dict, ["H"], coords, 'spherical')
     points = np.random.rand(10, 3)
 
     with pytest.raises(TypeError):
@@ -78,7 +78,7 @@ def test_evaluate_ehrenfest_force():
     """Test gbasis.evals.stress_tensor.evaluate_ehrenfest_force."""
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     coords = np.array([[0, 0, 0]])
-    basis = make_contractions(basis_dict, ["H"], coords)
+    basis = make_contractions(basis_dict, ["H"], coords, 'spherical')
     points = np.random.rand(10, 3)
 
     with pytest.raises(TypeError):
@@ -139,8 +139,8 @@ def test_evaluate_ehrenfest_hessian():
     """Test gbasis.evals.stress_tensor.evaluate_ehrenfest_hessian."""
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
     coords = np.array([[0, 0, 0]])
-    basis = make_contractions(basis_dict, ["H"], coords)
-    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps) for i in basis]
+    basis = make_contractions(basis_dict, ["H"], coords, 'spherical')
+    basis = [HortonContractions(i.angmom, i.coord, i.coeffs, i.exps, i.coord_type) for i in basis]
 
     points = np.random.rand(10, 3)
 
diff --git a/tests/test_wrappers.py b/tests/test_wrappers.py
index 70c6e077..89fbc38a 100644
--- a/tests/test_wrappers.py
+++ b/tests/test_wrappers.py
@@ -124,7 +124,7 @@ def test_from_pyscf():
     test = from_pyscf(mol)
 
     basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem"))
-    basis = make_contractions(basis_dict, ["Kr"], np.array([[1, 2, 3]]))
+    basis = make_contractions(basis_dict, ["Kr"], np.array([[1, 2, 3]]), 'spherical')
 
     with pytest.raises(ValueError):
 

From f7789135c2890d75bc2ba85ace85a72e95244d5f Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Thu, 11 Jan 2024 20:13:16 +0330
Subject: [PATCH 16/17] Remove support for mixed c/p coordinates from_pyscf

---
 gbasis/parsers.py  |  1 +
 gbasis/wrappers.py | 30 ++++++------------------------
 2 files changed, 7 insertions(+), 24 deletions(-)

diff --git a/gbasis/parsers.py b/gbasis/parsers.py
index 636670d3..3e63c6fb 100644
--- a/gbasis/parsers.py
+++ b/gbasis/parsers.py
@@ -210,6 +210,7 @@ def make_contractions(basis_dict, atoms, coords, coord_types):
 
     basis = []
     len_coord_types = len(coord_types)
+
     for atom, coord in zip(atoms, coords):
         for angmom, exps, coeffs in basis_dict[atom]:
             if type(coord_types) == str:
diff --git a/gbasis/wrappers.py b/gbasis/wrappers.py
index 50c776be..c94ec8ea 100644
--- a/gbasis/wrappers.py
+++ b/gbasis/wrappers.py
@@ -132,7 +132,9 @@ def angmom_components_sph(self):
 
         # pylint: disable=E1136
         basis.append(
-            IODataShell(angmom, mol.atcoords[shell.icenter], shell.coeffs, shell.exponents, shell.kinds[0])
+            IODataShell(
+                angmom, mol.atcoords[shell.icenter], shell.coeffs, shell.exponents, shell.kinds[0]
+            )
         )
 
     return basis
@@ -187,34 +189,14 @@ def angmom_components_sph(self):
             return super().angmom_components_sph
 
     basis = []
-    len_coord_types = len(coord_types)
     for atom, coord in mol._atom:
         basis_info = mol._basis[atom]
 
         for shell in basis_info:
             angmom = shell[0]
-
             exps_coeffs = np.vstack(shell[1:])
-            exps = exps_coeffs[:, 0]
-
-            coeffs = exps_coeffs[:, 1:]
-
-            if type(coord_types) == str:
-            # if coord_types given as a single string, assign the specified type to all contractions for all atoms
-                if coord_types == "spherical":
-                    basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps), 'p'))
-                elif coord_types == "cartesian":
-                    basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps), 'c'))
-                else:
-                    raise ValueError("If coord_types is a string, it must be either 'spherical' or 'cartesian'.")
-            elif type(coord_types) == list:
-                # if coord_types given as a list, assign the specified type to each atom's contractions individually
-                if len_coord_types == mol.nbas:
-                    basis.append(PyscfShell(angmom, np.array(coord), np.array(coeffs), np.array(exps), coord_types.pop(0)))
-                else:
-                    raise ValueError(
-                        "If coord_types is a list, it must be the same length as the total number of contractions.")
-            else:
-                raise TypeError("coord_types must be a string or list of strings.")
+            exps = np.array(exps_coeffs[:, 0])
+            coeffs = np.array(exps_coeffs[:, 1:])
+            basis.append(PyscfShell(angmom, np.array(coord), coeffs, exps, coord_types))
 
     return tuple(basis)

From 449c21f24091e23d780804e4b28a15a06fa37821 Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Thu, 11 Jan 2024 20:30:24 +0330
Subject: [PATCH 17/17] Refactor make_contractions

---
 gbasis/parsers.py | 39 +++++++++++++++++++++------------------
 1 file changed, 21 insertions(+), 18 deletions(-)

diff --git a/gbasis/parsers.py b/gbasis/parsers.py
index 3e63c6fb..3b28189f 100644
--- a/gbasis/parsers.py
+++ b/gbasis/parsers.py
@@ -209,25 +209,28 @@ def make_contractions(basis_dict, atoms, coords, coord_types):
         raise ValueError("Number of atoms must be equal to the number of rows in the coordinates.")
 
     basis = []
-    len_coord_types = len(coord_types)
+    # expected number of coordinates
+    num_coord_types = sum([len(basis_dict[i]) for i in atoms])
+
+    # check and assign coord_types
+    if isinstance(coord_types, str):
+        if coord_types not in ["c", "cartesian", "p", "spherical"]:
+            raise ValueError(
+                f"If coord_types is a string, it must be either 'spherical'/'p' or 'cartesian'/'c'."
+                f"got {coord_types}"
+            )
+        coord_types = [coord_types] * num_coord_types
+
+    if len(coord_types) != num_coord_types:
+        raise ValueError(
+            f"If coord_types is a list, it must be the same length as the total number of contractions."
+            f"got {len(coord_types)}"
+        )
 
+    # make shells
     for atom, coord in zip(atoms, coords):
         for angmom, exps, coeffs in basis_dict[atom]:
-            if type(coord_types) == str:
-                # if coord_types given as a single string, assign the specified type to all contractions for all atoms
-                if coord_types == "spherical" or coord_types == "p":
-                    basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, 'p'))
-                elif coord_types == "cartesian" or coord_types == "c":
-                    basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, 'c'))
-                else:
-                    raise ValueError("If coord_types is a string, it must be either 'spherical'/'p' or 'cartesian'/'c'.")
-            elif type(coord_types) == list:
-                # if coord_types given as a list, assign the specified type to each atom's contractions individually
-                if len_coord_types == sum([len(basis_dict[i]) for i in atoms]):
-                    basis.append(GeneralizedContractionShell(angmom, coord, coeffs, exps, coord_types.pop(0)))
-                else:
-                    raise ValueError("If coord_types is a list, it must be the same length as the total number of contractions.")
-            else:
-                raise TypeError("coord_types must be a string or list of strings.")
-
+            basis.append(
+                GeneralizedContractionShell(angmom, coord, coeffs, exps, coord_types.pop(0))
+            )
     return tuple(basis)