diff --git a/.editorconfig b/.editorconfig
index df64623b..dde5c5bf 100644
--- a/.editorconfig
+++ b/.editorconfig
@@ -1,4 +1,4 @@
-# EditorConfig is awesome: https://EditorConfig.org
+# EditorConfig: https://EditorConfig.org
 
 root = true
 
diff --git a/gbasis/base_four_symm.py b/gbasis/base_four_symm.py
index 28f375a6..7ce64df8 100644
--- a/gbasis/base_four_symm.py
+++ b/gbasis/base_four_symm.py
@@ -2,9 +2,10 @@
 import abc
 import itertools as it
 
+import numpy as np
+
 from gbasis.base import BaseGaussianRelatedArray
 from gbasis.spherical import generate_transformation
-import numpy as np
 
 
 # pylint: disable=W0235
diff --git a/gbasis/base_one.py b/gbasis/base_one.py
index 5da708c3..1d53d802 100644
--- a/gbasis/base_one.py
+++ b/gbasis/base_one.py
@@ -1,9 +1,10 @@
 """Base class for arrays that depend on one contracted Gaussian."""
 import abc
 
+import numpy as np
+
 from gbasis.base import BaseGaussianRelatedArray
 from gbasis.spherical import generate_transformation
-import numpy as np
 
 
 # pylint: disable=W0235
diff --git a/gbasis/base_two_asymm.py b/gbasis/base_two_asymm.py
index 5d5c2e69..edcc709c 100644
--- a/gbasis/base_two_asymm.py
+++ b/gbasis/base_two_asymm.py
@@ -1,9 +1,10 @@
 """Base class for arrays that depend on two contracted Gaussians."""
 import abc
 
+import numpy as np
+
 from gbasis.base import BaseGaussianRelatedArray
 from gbasis.spherical import generate_transformation
-import numpy as np
 
 
 # pylint: disable=W0235
diff --git a/gbasis/base_two_symm.py b/gbasis/base_two_symm.py
index 40e0cdeb..7dd3f80b 100644
--- a/gbasis/base_two_symm.py
+++ b/gbasis/base_two_symm.py
@@ -2,9 +2,10 @@
 
 import abc
 
+import numpy as np
+
 from gbasis.base import BaseGaussianRelatedArray
 from gbasis.spherical import generate_transformation
-import numpy as np
 
 
 # pylint: disable=W0235
diff --git a/gbasis/contractions.py b/gbasis/contractions.py
index 1d4ee789..21d5af4a 100644
--- a/gbasis/contractions.py
+++ b/gbasis/contractions.py
@@ -1,7 +1,9 @@
 """Data class for contractions of Gaussian-type primitives."""
 
 from numbers import Integral
+
 import numpy as np
+
 from gbasis.utils import factorial2
 
 
@@ -31,7 +33,8 @@ class GeneralizedContractionShell:
     .. math::
 
         \phi (\mathbf{r} | \mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) =
-        N_{\phi} (\mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) \sum_i d_i g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a})
+        N_{\phi} (\mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha})
+        \sum_i d_i g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a})
 
     where :math:`d_i` is the contraction coefficient of the primitive and :math:`N_{\phi}` is the
     normalization constant of the contraction.
@@ -63,7 +66,8 @@ class GeneralizedContractionShell:
     .. math::
 
         \phi_j (\mathbf{r} | \mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) =
-        N_{\phi} (\mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) \sum_i d_{ij} g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a})
+        N_{\phi} (\mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha})
+        \sum_i d_{ij} g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a})
 
     Attributes
     ----------
@@ -427,8 +431,8 @@ def angmom_components_sph(self):
             return tuple(["c1", "s1", "c0"])
         else:
             return tuple(
-                ["s{}".format(m) for m in range(self.angmom, 0, -1)]
-                + ["c{}".format(m) for m in range(self.angmom + 1)]
+                [f"s{m}" for m in range(self.angmom, 0, -1)]
+                + [f"c{m}" for m in range(self.angmom + 1)]
             )
 
     @property
diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py
index 82341804..214512d9 100644
--- a/gbasis/evals/density.py
+++ b/gbasis/evals/density.py
@@ -1,9 +1,10 @@
 """Density Evaluation."""
-from gbasis.evals.eval import evaluate_basis
-from gbasis.evals.eval_deriv import evaluate_deriv_basis
 import numpy as np
 from scipy.special import comb
 
+from gbasis.evals.eval import evaluate_basis
+from gbasis.evals.eval_deriv import evaluate_deriv_basis
+
 
 def evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval):
     """Return the evaluation of the density given the evaluated orbitals.
@@ -218,15 +219,16 @@ def evaluate_deriv_density(
             \sum_{l_x=0}^{L_x} \sum_{l_y=0}^{L_y} \sum_{l_z=0}^{L_z}
             \binom{L_x}{l_x} \binom{L_y}{l_y} \binom{L_z}{l_z}
             \sum_{ij} \gamma_{ij}
-            \frac{\partial^{l_x + l_y + l_z} \rho(\mathbf{r})}{\partial x^{l_x} \partial y^{l_y} \partial z^{l_z}}
+            \frac{\partial^{l_x + l_y + l_z}
+            \rho(\mathbf{r})}{\partial x^{l_x} \partial y^{l_y} \partial z^{l_z}}
             \frac{
                 \partial^{L_x + L_y + L_z - l_x - l_y - l_z} \rho(\mathbf{r})
             }{
                 \partial x^{L_x - l_x} \partial y^{L_y - l_y} \partial z^{L_z - l_z}
             }
 
-    where :math:`L_x, L_y, L_z` are the orders of the derivative relative to the :math:`x, y, \text{and} z` components,
-    respectively.
+    where :math:`L_x, L_y, L_z` are the orders of the derivative relative to the
+    :math:`x, y, \text{and} z` components, respectively.
 
     Parameters
     ----------
@@ -320,7 +322,7 @@ def evaluate_density_gradient(
             \frac{\partial}{\partial x} \rho(\mathbf{r})\\\\
             \frac{\partial}{\partial y} \rho(\mathbf{r})\\\\
             \frac{\partial}{\partial z} \rho(\mathbf{r})
-            \end{bmatrix}      
+            \end{bmatrix}
 
     Parameters
     ----------
@@ -489,7 +491,7 @@ def evaluate_density_hessian(
     r"""Return the Hessian of the density evaluated at the given points.
 
     .. math::
-        
+
         H[\rho(\mathbf{r})]
         =
         \begin{bmatrix}
diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py
index 4d850975..a13cc9d8 100644
--- a/gbasis/evals/electrostatic_potential.py
+++ b/gbasis/evals/electrostatic_potential.py
@@ -1,7 +1,8 @@
 """Module for computing electrostatic potential integrals."""
-from gbasis.integrals.point_charge import point_charge_integral
 import numpy as np
 
+from gbasis.integrals.point_charge import point_charge_integral
+
 
 def electrostatic_potential(
     basis,
@@ -18,12 +19,13 @@ def electrostatic_potential(
 
         - \left(
             - \sum_A \frac{Z_A}{|\mathbf{R}_C - \mathbf{R}_A|}
-            + \sum_{ab} \gamma_{ab} \int \phi_a(\mathbf{r}) \frac{-1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+            + \sum_{ab} \gamma_{ab} \int \phi_a(\mathbf{r}) \frac{-1}{|\mathbf{r}
+            - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
         \right)
 
-    where :math:`\mathbf{R}_C` is the coordinate of a unitary point charge, :math:`\mathbf{R}_A` is the
-    coordinate of the nucleus :math:`A`, :math:`Z_A` its charge, and :math:`\gamma_{ab}` is the
-    one-electron density matrix.
+    where :math:`\mathbf{R}_C` is the coordinate of a unitary point charge, :math:`\mathbf{R}_A`
+    is the coordinate of the nucleus :math:`A`, :math:`Z_A` its charge, and :math:`\gamma_{ab}`
+    is the one-electron density matrix.
 
     Parameters
     ----------
diff --git a/gbasis/evals/eval.py b/gbasis/evals/eval.py
index bd14b8d6..73f29952 100644
--- a/gbasis/evals/eval.py
+++ b/gbasis/evals/eval.py
@@ -1,8 +1,9 @@
 """Functions for evaluating Gaussian contractions."""
+import numpy as np
+
 from gbasis.base_one import BaseOneIndex
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.evals._deriv import _eval_deriv_contractions
-import numpy as np
 
 
 class Eval(BaseOneIndex):
diff --git a/gbasis/evals/eval_deriv.py b/gbasis/evals/eval_deriv.py
index 245aa1cb..3ee85c2d 100644
--- a/gbasis/evals/eval_deriv.py
+++ b/gbasis/evals/eval_deriv.py
@@ -1,9 +1,12 @@
 """Functions for evaluating Gaussian primitives."""
+import numpy as np
+
 from gbasis.base_one import BaseOneIndex
 from gbasis.contractions import GeneralizedContractionShell
-from gbasis.evals._deriv import _eval_deriv_contractions
-from gbasis.evals._deriv import _eval_first_second_order_deriv_contractions
-import numpy as np
+from gbasis.evals._deriv import (
+    _eval_deriv_contractions,
+    _eval_first_second_order_deriv_contractions,
+)
 
 
 class EvalDeriv(BaseOneIndex):
@@ -150,13 +153,13 @@ def evaluate_deriv_basis(
     The derivative (to arbitrary orders) of a basis function is given by:
 
     .. math::
-    
+
         \frac{\partial^{m_x + m_y + m_z}}{\partial x^{m_x} \partial y^{m_y} \partial z^{m_z}}
         \phi (\mathbf{r})
-    
-    where :math:`m_x, m_y, m_z` are the orders of the derivative with respect to x, y, and z, 
-    :math:`\phi` is the basis function (a generalized contraction shell), and :math:`\mathbf{r}_n` are
-    the coordinate of the points in space where the basis function is evaluated.
+
+    where :math:`m_x, m_y, m_z` are the orders of the derivative with respect to x, y, and z,
+    :math:`\phi` is the basis function (a generalized contraction shell), and :math:`\mathbf{r}_n`
+    are the coordinate of the points in space where the basis function is evaluated.
 
     Parameters
     ----------
diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py
index 6a33ac51..b94e0b91 100644
--- a/gbasis/evals/stress_tensor.py
+++ b/gbasis/evals/stress_tensor.py
@@ -1,10 +1,11 @@
 """Module for computing properties related to the stress tensor."""
+import numpy as np
+
 from gbasis.evals.density import (
     evaluate_density_laplacian,
     evaluate_deriv_density,
     evaluate_deriv_reduced_density_matrix,
 )
-import numpy as np
 
 
 # TODO: need to be tested against reference
@@ -132,31 +133,31 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0,
 
     .. math::
 
-            F_{j}(\mathbf{r} | \alpha, \beta)
-            =&- \sum_i \frac{\partial}{\partial r_i} \boldsymbol{\sigma}_{ij}\\
-            =&
-            \alpha
-            \sum_i
-            \left.
-              \frac{\partial^3}{\partial r^2_i \partial r'_j} \gamma(\mathbf{r}, \mathbf{r}')
-            \right|_{\mathbf{r} = \mathbf{r}'}\\
-            &- (1 - \alpha)
-            \sum_i
-            \left.
-              \frac{\partial^3}{\partial r^2_i \partial r_j} \gamma(\mathbf{r}, \mathbf{r})
-            \right|_{\mathbf{r} = \mathbf{r}'}
-            - (1 - 2\alpha)
-            \sum_i
-            \left.
-              \frac{\partial^3}{\partial r_i \partial r_j \partial r'_i} \gamma(\mathbf{r}, \mathbf{r})
-            \right|_{\mathbf{r} = \mathbf{r}'}\\
-            &+ \frac{1}{2} \beta
-            \left(
-            \frac{\partial^3}{\partial r_j \partial x^2}
-            + \frac{\partial^3}{\partial r_j \partial y^2}
-            + \frac{\partial^3}{\partial r_j \partial z^2}
-            \right)
-            \rho(\mathbf{r})
+        F_{j}(\mathbf{r} | \alpha, \beta)
+        =&- \sum_i \frac{\partial}{\partial r_i} \boldsymbol{\sigma}_{ij}\\
+        =&
+        \alpha
+        \sum_i
+        \left.
+          \frac{\partial^3}{\partial r^2_i \partial r'_j} \gamma(\mathbf{r}, \mathbf{r}')
+        \right|_{\mathbf{r} = \mathbf{r}'}\\
+        &- (1 - \alpha)
+        \sum_i
+        \left.
+          \frac{\partial^3}{\partial r^2_i \partial r_j} \gamma(\mathbf{r}, \mathbf{r})
+        \right|_{\mathbf{r} = \mathbf{r}'}
+        - (1 - 2\alpha)
+        \sum_i
+        \left.
+          \frac{\partial^3}{\partial r_i \partial r_j \partial r'_i} \gamma(\mathbf{r}, \mathbf{r})
+        \right|_{\mathbf{r} = \mathbf{r}'}\\
+        &+ \frac{1}{2} \beta
+        \left(
+        \frac{\partial^3}{\partial r_j \partial x^2}
+        + \frac{\partial^3}{\partial r_j \partial y^2}
+        + \frac{\partial^3}{\partial r_j \partial z^2}
+        \right)
+        \rho(\mathbf{r})
 
     Parameters
     ----------
diff --git a/gbasis/integrals/_diff_operator_int.py b/gbasis/integrals/_diff_operator_int.py
index d590d8bf..a01be2f5 100644
--- a/gbasis/integrals/_diff_operator_int.py
+++ b/gbasis/integrals/_diff_operator_int.py
@@ -1,9 +1,10 @@
 """Integrals over differential operator involving contracted Cartesian Gaussians."""
+import numpy as np
+
 from gbasis.integrals._moment_int import (
     _cleanup_intermediate_integrals,
     _compute_multipole_moment_integrals_intermediate,
 )
-import numpy as np
 
 
 # FIXME: returns nan when exponent is zero
diff --git a/gbasis/integrals/_one_elec_int.py b/gbasis/integrals/_one_elec_int.py
index 688d7d9d..a3691be4 100644
--- a/gbasis/integrals/_one_elec_int.py
+++ b/gbasis/integrals/_one_elec_int.py
@@ -1,5 +1,6 @@
 """One-electron integrals involving Contracted Cartesian Gaussians."""
 import numpy as np
+
 from gbasis.utils import factorial2
 
 
diff --git a/gbasis/integrals/_two_elec_int.py b/gbasis/integrals/_two_elec_int.py
index 6ff2a651..526e44d5 100644
--- a/gbasis/integrals/_two_elec_int.py
+++ b/gbasis/integrals/_two_elec_int.py
@@ -1,5 +1,6 @@
 """Two-electron integrals involving Contracted Cartesian Gaussians."""
 import numpy as np
+
 from gbasis.utils import factorial2
 
 # pylint: disable=C0103,R0914,R0915
diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py
index e31986ef..74275e65 100644
--- a/gbasis/integrals/angular_momentum.py
+++ b/gbasis/integrals/angular_momentum.py
@@ -1,11 +1,12 @@
 """Module for evaluating the integral over the angular momentum operator."""
+import numpy as np
+
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._diff_operator_int import (
     _compute_differential_operator_integrals_intermediate,
 )
 from gbasis.integrals._moment_int import _compute_multipole_moment_integrals_intermediate
-import numpy as np
 
 
 # TODO: need to test against reference
@@ -161,17 +162,18 @@ def angular_momentum_integral(basis, transform=None):
 
     .. math::
 
-            \left< \hat{\mathbf{L}} \right>
-            &= \int \phi_a(\mathbf{r}) \left( -i \mathbf{r} \times \nabla \right) \phi_b(\mathbf{r}) d\mathbf{r}\\
-            &= -i
-            \begin{bmatrix}
-            \int \phi_a(\mathbf{r}) y\frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}
-            - \int \phi_a(\mathbf{r}) z\frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
-            \int \phi_a(\mathbf{r}) z\frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}
-            - \int \phi_a(\mathbf{r}) x\frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
-            \int \phi_a(\mathbf{r}) x\frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}
-            - \int \phi_a(\mathbf{r}) y\frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
-            \end{bmatrix}
+        \left< \hat{\mathbf{L}} \right>
+        &= \int \phi_a(\mathbf{r}) \left( -i \mathbf{r} \times \nabla \right) \phi_b(\mathbf{r})
+        d\mathbf{r}\\
+        &= -i
+        \begin{bmatrix}
+        \int \phi_a(\mathbf{r}) y\frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}
+        - \int \phi_a(\mathbf{r}) z\frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+        \int \phi_a(\mathbf{r}) z\frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}
+        - \int \phi_a(\mathbf{r}) x\frac{\partial}{\partial z} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+        \int \phi_a(\mathbf{r}) x\frac{\partial}{\partial y} \phi_b(\mathbf{r}) d\mathbf{r}
+        - \int \phi_a(\mathbf{r}) y\frac{\partial}{\partial x} \phi_b(\mathbf{r}) d\mathbf{r}\\\\
+        \end{bmatrix}
 
     Parameters
     ----------
diff --git a/gbasis/integrals/electron_repulsion.py b/gbasis/integrals/electron_repulsion.py
index 4643bbbe..de824e61 100644
--- a/gbasis/integrals/electron_repulsion.py
+++ b/gbasis/integrals/electron_repulsion.py
@@ -1,4 +1,6 @@
 """Electron-electron repulsion integral."""
+import numpy as np
+
 from gbasis.base_four_symm import BaseFourIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._two_elec_int import (
@@ -6,7 +8,6 @@
     _compute_two_elec_integrals_angmom_zero,
 )
 from gbasis.integrals.point_charge import PointChargeIntegral
-import numpy as np
 
 
 class ElectronRepulsionIntegral(BaseFourIndexSymmetric):
diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py
index d944bd47..f46511b1 100644
--- a/gbasis/integrals/kinetic_energy.py
+++ b/gbasis/integrals/kinetic_energy.py
@@ -1,8 +1,9 @@
 """Module for evaluating the kinetic energy integral."""
+import numpy as np
+
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals
-import numpy as np
 
 
 class KineticEnergyIntegral(BaseTwoIndexSymmetric):
@@ -126,7 +127,8 @@ def kinetic_energy_integral(basis, transform=None):
 
         \begin{split}
         \left< \hat{T} \right>
-        &= \int \phi_a(\mathbf{r}) \left( -\frac{1}{2} \nabla^2 \right) \phi_b(\mathbf{r}) d\mathbf{r}\\
+        &= \int \phi_a(\mathbf{r}) \left( -\frac{1}{2} \nabla^2 \right) \phi_b(\mathbf{r})
+        d\mathbf{r}\\
         &= -\frac{1}{2}
         \left(
         \int \phi_a(\mathbf{r}) \frac{\partial^2}{\partial x^2} \phi_b(\mathbf{r}) d\mathbf{r}
diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py
index a6df60ff..a16f4ae4 100644
--- a/gbasis/integrals/libcint.py
+++ b/gbasis/integrals/libcint.py
@@ -2,26 +2,15 @@
 Python C-API bindings for ``libcint`` GTO integrals library.
 
 """
-
+import re
 from contextlib import contextmanager
-
-from ctypes import CDLL, POINTER, Structure, cdll, byref, c_int, c_double, c_void_p
-
-from itertools import chain
-
+from ctypes import CDLL, POINTER, Structure, byref, c_double, c_int, cdll
 from operator import attrgetter
-
 from pathlib import Path
 
-import re
-
 import numpy as np
-
 from scipy.special import factorial
 
-from gbasis.utils import factorial2
-
-
 __all__ = [
     "LIBCINT",
     "CBasis",
@@ -218,7 +207,7 @@ class _LibCInt:
 
     """
 
-    _libcint: CDLL = cdll.LoadLibrary((Path(__file__).parent / "lib" / "libcint.so"))
+    _libcint: CDLL = cdll.LoadLibrary(Path(__file__).parent / "lib" / "libcint.so")
     r"""
     ``libcint`` shared object library.
 
@@ -235,7 +224,7 @@ def __new__(cls):
 
         """
         if not hasattr(cls, "_instance"):
-            cls._instance = super(_LibCInt, cls).__new__(cls)
+            cls._instance = super().__new__(cls)
         return cls._instance
 
     def __init__(self):
@@ -386,9 +375,11 @@ class CBasis:
 
     Methods
     -------
-    make_int1e(self, func_name, components=tuple(), constant=None, is_complex=False, origin=False, inv_origin=False)
+    make_int1e(self, func_name, components=tuple(), constant=None,
+               is_complex=False, origin=False, inv_origin=False)
         Make an instance-bound 1-electron integral method from a ``libcint`` function.
-    make_int2e(self, func_name, components=tuple(), constant=None, is_complex=False, origin=False, inv_origin=False)
+    make_int2e(self, func_name, components=tuple(), constant=None,
+               is_complex=False, origin=False, inv_origin=False)
         Make an instance-bound 2-electron integral method from a ``libcint`` function.
     overlap(self)
         Compute the overlap integrals.
@@ -630,12 +621,13 @@ def make_int1e(
         components : tuple, default=()
             Shape of components in each integral element.
             E.g., for normal integrals, ``comp=()``, while for nuclear gradients,
-            ``components=(Natm, 3)``, and for nuclear Hessians, ``components=(Natm, Natm, 3, 3)``, etc.
+            ``components=(Natm, 3)``, and for nuclear Hessians,
+            ``components=(Natm, Natm, 3, 3)``, etc.
         constant : (float | complex), default=1.
             A constant by which to multiply the whole integral array.
         is_complex : bool, default=False
             Whether the components in each integral element are complex. Not required if only
-            multiplying by a complex constant using the ``constant`` keyword argument..
+            multiplying by a complex constant using the ``constant`` keyword argument.
         origin : bool, default=False
             Whether you must specify an origin ``R`` for the integral computation.
         inv_origin : bool, default=False
@@ -656,7 +648,7 @@ def make_int1e(
         if is_complex:
             components += (2,)
         prod_comp = np.prod(components, dtype=int)
-        out_shape = (self.nbfn, self.nbfn) + components
+        out_shape = (self.nbfn, self.nbfn, *components)
         buf_shape = prod_comp * self._max_off**2
 
         # Handle [inv_]origin argument (prevent shadowing)
@@ -789,7 +781,8 @@ def make_int2e(
         components : tuple, default=()
             Shape of components in each integral element.
             E.g., for normal integrals, ``components=(1,)``, while for nuclear gradients,
-            ``components=(Natm, 3)``, and for nuclear Hessians, ``components=(Natm, Natm, 3, 3)``, etc.
+            ``components=(Natm, 3)``, and for nuclear Hessians,
+            ``components=(Natm, Natm, 3, 3)``, etc.
         constant : (float | complex), default=1.
             A constant by which to multiply the whole integral array.
         is_complex : bool, default=False
@@ -1001,8 +994,8 @@ def overlap_integral(self, notation="physicist", transform=None):
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1024,8 +1017,8 @@ def kinetic_energy_integral(self, notation="physicist", transform=None):
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1047,8 +1040,8 @@ def nuclear_attraction_integral(self, notation="physicist", transform=None):
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1070,8 +1063,8 @@ def electron_repulsion_integral(self, notation="physicist", transform=None):
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1095,8 +1088,8 @@ def r_inv_integral(self, origin=None, notation="physicist", transform=None):
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1120,8 +1113,8 @@ def momentum_integral(self, origin=None, notation="physicist", transform=None):
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1145,8 +1138,8 @@ def angular_momentum_integral(self, origin=None, notation="physicist", transform
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1175,8 +1168,8 @@ def point_charge_integral(
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
@@ -1210,8 +1203,8 @@ def moment_integral(self, orders, origin=None, notation="physicist", transform=N
         notation : ("physicist" | "chemist"), default="physicist"
             Axis order convention.
         transform : np.ndarray(K, K_cont)
-            Transformation matrix from the basis set in the given coordinate system (e.g. AO) to linear
-            combinations of contractions (e.g. MO).
+            Transformation matrix from the basis set in the given coordinate system (e.g. AO)
+            to linear combinations of contractions (e.g. MO).
             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.
diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py
index 78e486a7..7b04c8f7 100644
--- a/gbasis/integrals/moment.py
+++ b/gbasis/integrals/moment.py
@@ -1,8 +1,9 @@
 """Module for computing the moments of a basis set."""
+import numpy as np
+
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._moment_int import _compute_multipole_moment_integrals
-import numpy as np
 
 
 class Moment(BaseTwoIndexSymmetric):
@@ -162,10 +163,11 @@ def moment_integral(basis, moment_coord, moment_orders, transform=None):
 
     .. math::
 
-        \int \phi_a (\mathbf{r}) (x - X_C)^{c_x} (y - Y_C)^{c_y} (z - Z_C)^{c_z} \phi_b (\mathbf{r}) d\mathbf{r}
+        \int \phi_a (\mathbf{r}) (x - X_C)^{c_x} (y - Y_C)^{c_y} (z - Z_C)^{c_z}
+        \phi_b (\mathbf{r}) d\mathbf{r}
 
-    where :math:`X_C`, :math:`Y_C`, and :math:`Z_C` are the coordinates of the center of the moment, and
-    :math:`c_x`, :math:`c_y`, and :math:`c_z` are the orders of the moment.
+    where :math:`X_C`, :math:`Y_C`, and :math:`Z_C` are the coordinates of the center of the moment,
+    and :math:`c_x`, :math:`c_y`, and :math:`c_z` are the orders of the moment.
 
     Parameters
     ----------
diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py
index 57d3f4eb..93ff0fa2 100644
--- a/gbasis/integrals/momentum.py
+++ b/gbasis/integrals/momentum.py
@@ -1,8 +1,9 @@
 """Module for evaluating the integral over the momentum operator."""
+import numpy as np
+
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals
-import numpy as np
 
 
 # TODO: need to test against reference
diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py
index 00e2241e..b143e16c 100644
--- a/gbasis/integrals/nuclear_electron_attraction.py
+++ b/gbasis/integrals/nuclear_electron_attraction.py
@@ -1,7 +1,8 @@
 """Module for computing the nuclear electron attraction."""
-from gbasis.integrals.point_charge import point_charge_integral
 import numpy as np
 
+from gbasis.integrals.point_charge import point_charge_integral
+
 
 def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, transform=None):
     """Return the nuclear electron attraction integrals of the basis set in the Cartesian form.
@@ -11,7 +12,7 @@ def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges,
         \sum_C \int \phi_a(\mathbf{r}) \frac{-Z_C}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r})
     d\mathbf{r}
 
-    where :math:`\mathbf{R}_C` is the coordinate of the nucleus :math:`C` and :math:`Z_C` its charge.
+    where :math:`\mathbf{R}_C` is the coordinate of the nucleus :math:`C` and :math:`Z_C` its charge
 
     Parameters
     ----------
diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py
index 068959d0..54d23523 100644
--- a/gbasis/integrals/overlap.py
+++ b/gbasis/integrals/overlap.py
@@ -1,9 +1,9 @@
 """Functions for computing overlap of a basis set."""
+import numpy as np
 
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._moment_int import _compute_multipole_moment_integrals
-import numpy as np
 
 
 class Overlap(BaseTwoIndexSymmetric):
@@ -62,10 +62,10 @@ def construct_array_contraction(contractions_one, contractions_two, tol_screen=N
             Contracted Cartesian Gaussians (of the same shell) associated with the second index of
             the overlap.
         tol_screen : None or float, optional
-            The tolerance used for screening overlap integrals. The `tol_screen` is combined with the
-            minimum contraction exponents to compute a cutoff which is compared to the distance between
-            the contraction centers to decide whether the overlap integral should be set to zero (i.e.,
-            screened). If `None`, no screening is performed.
+            The tolerance used for screening overlap integrals. The `tol_screen` is combined with
+            the minimum contraction exponents to compute a cutoff which is compared to the distance
+            between the contraction centers to decide whether the overlap integral should be set to
+            zero (i.e.,screened). If `None`, no screening is performed.
 
         Returns
         -------
@@ -181,7 +181,8 @@ def is_integral_screened(contractions_one, contractions_two, tol_screen):
     where :math:`d` is the cut-off distance at which shells do not interact with each other.
     :math:`A` and `B` are the atoms each contraction are respectively centered on.
     :math: `\alpha` is the gaussian exponent
-    :math: `s` and `t` index the primitive Gaussian shells centered on atom `A` and `B` respectively.
+    :math: `s` and `t` index the primitive Gaussian shells centered on atom `A` and `B`
+    respectively.
 
     Parameters
     ----------
diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py
index 294f1bd2..0dd1cab6 100644
--- a/gbasis/integrals/point_charge.py
+++ b/gbasis/integrals/point_charge.py
@@ -1,9 +1,10 @@
 """Module for computing point charge integrals."""
+import numpy as np
+from scipy.special import hyp1f1  # pylint: disable=E0611
+
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._one_elec_int import _compute_one_elec_integrals
-import numpy as np
-from scipy.special import hyp1f1  # pylint: disable=E0611
 
 
 class PointChargeIntegral(BaseTwoIndexSymmetric):
@@ -37,7 +38,8 @@ class PointChargeIntegral(BaseTwoIndexSymmetric):
         Boys function used to evaluate the one-electron integral.
         `M` is the number of orders that will be evaluated.
         `K_a` and `K_b` are the number of primitives on the left and right side, respectively.
-    construct_array_contraction(self, contractions_one, contractions_two, points_coords, points_charge)
+    construct_array_contraction(self, contractions_one, contractions_two,
+                                points_coords, points_charge)
         Return the point charge integrals for the given `GeneralizedContractionShell` instances.
         `M_1` is the number of segmented contractions with the same exponents (and angular
         momentum) associated with the first index.
@@ -271,11 +273,12 @@ def point_charge_integral(basis, points_coords, points_charge, transform=None):
 
     .. math::
 
-        V_{ab} = \int \phi_a(\mathbf{r}) \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+        V_{ab} = \int \phi_a(\mathbf{r})
+        \frac{1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
 
     where :math:`\mathbf{R}_C` is the position of the point charge :math:`C` and :math:`V_{ab}` is
-    the interaction integral between the pair of basis functions :math:`\phi_a` and :math:`\phi_b`, and
-    the point charge.
+    the interaction integral between the pair of basis functions :math:`\phi_a` and :math:`\phi_b`,
+    and the point charge.
 
     Parameters
     ----------
diff --git a/gbasis/parsers.py b/gbasis/parsers.py
index 686412f6..5e3f5ab2 100644
--- a/gbasis/parsers.py
+++ b/gbasis/parsers.py
@@ -2,9 +2,10 @@
 
 import re
 
-from gbasis.contractions import GeneralizedContractionShell
 import numpy as np
 
+from gbasis.contractions import GeneralizedContractionShell
+
 
 def parse_nwchem(nwchem_basis_file):
     """Parse nwchem basis set file.
@@ -27,7 +28,7 @@ def parse_nwchem(nwchem_basis_file):
 
     """
     # pylint: disable=R0914
-    with open(nwchem_basis_file, "r") as basis_fh:
+    with open(nwchem_basis_file) as basis_fh:
         nwchem_basis = basis_fh.read()
 
     data = re.split(r"\n\s*(\w[\w]?)[ ]+(\w+)\s*\n", nwchem_basis)
@@ -177,13 +178,14 @@ def make_contractions(basis_dict, atoms, coords, coord_types):
         Atoms at which the contractions are centered.
     coords : np.ndarray(N, 3)
         Coordinates of each atom.
-    coord_types : {"cartesian"/"c", list/tuple of "cartesian"/"c" or "spherical"/"p", "spherical"/"p"}
+    coord_types : {"cartesian"/"c", list/tuple of "cartesian"/"c" or
+                   "spherical"/"p", "spherical"/"p"}
         Types of the coordinate system for the contractions.
         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".
+        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".
 
     Returns
     -------
@@ -227,8 +229,8 @@ def make_contractions(basis_dict, atoms, coords, 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)}"
+            f"If coord_types is a list, it must be the same length as the total number "
+            f"of contractions. got {len(coord_types)}"
         )
 
     # make shells
diff --git a/gbasis/spherical.py b/gbasis/spherical.py
index 610ef38a..7137e2a0 100644
--- a/gbasis/spherical.py
+++ b/gbasis/spherical.py
@@ -4,6 +4,7 @@
 
 import numpy as np
 from scipy.special import comb, factorial
+
 from gbasis.utils import factorial2
 
 
@@ -45,7 +46,8 @@ def expansion_coeff(angmom, mag, i, j, k):
         C^{angmom,mag,i,j,k} = -1^{i + k - \text{shift_factor}} * (1/4)^i * {angmom \choose i}
         * {(angmom - i) \choose (|mag| + i)} * {i \choose j} * {|mag| \choose 2 * k},
 
-    where :math:`shift_factor = 0` if :math:`mag >= 0` and :math:`shift_factor = 1/2` if :math:`mag < 0`    .
+    where :math:`shift_factor = 0` if :math:`mag >= 0` and
+    :math:`shift_factor = 1/2` if :math:`mag < 0`.
 
     Parameters
     ----------
@@ -122,7 +124,8 @@ def harmonic_norm(angmom, mag):
     Calculate the normalization constant of a real solid harmonic.
 
     .. math::
-        N^S_{angmom,m} = 1/(2^{|m|} * angmom!) * \sqrt{2 * (angmom + |{m}|)! * (angmom - |{m}|)! / 2^{\delta_0^m}},
+        N^S_{angmom,m} = 1/(2^{|m|} * angmom!) * \sqrt{2 * (angmom + |{m}|)! * (angmom - |{m}|)! /
+        2^{\delta_0^m}},
 
     where :math: `del{0, m}` is the Kronecker delta of 0 and m.
 
@@ -300,10 +303,7 @@ def generate_transformation(angmom, cartesian_order, spherical_order, apply_from
         len(spherical_order) == 2 * angmom + 1
         # Strip out "-" from the ordering to make sure the right components are there
         and {x.replace("-", "") for x in spherical_order}
-        == set(
-            ["s{}".format(m) for m in range(angmom, 0, -1)]
-            + ["c{}".format(m) for m in range(angmom + 1)]
-        )
+        == set([f"s{m}" for m in range(angmom, 0, -1)] + [f"c{m}" for m in range(angmom + 1)])
     ):
         raise ValueError(
             "`spherical_order` must contain exactly 2 * `angmom` + 1 pure function strings, with "
diff --git a/gbasis/wrappers.py b/gbasis/wrappers.py
index cb550b2f..6cc8e380 100644
--- a/gbasis/wrappers.py
+++ b/gbasis/wrappers.py
@@ -1,8 +1,7 @@
 """Module for interfacing to other quantum chemistry packages."""
-
-from gbasis.contractions import GeneralizedContractionShell
 import numpy as np
 
+from gbasis.contractions import GeneralizedContractionShell
 
 CONVENTIONS_LIBCINT = {
     (5, "p"): ["s5", "s4", "s3", "s2", "s1", "c0", "c1", "c2", "c3", "c4", "c5"],
@@ -186,8 +185,8 @@ def angmom_components_sph(self):
             """
             if self.angmom not in sph_conventions:
                 raise ValueError(
-                    "Given convention does not support spherical contractions for the angular "
-                    "momentum {0}".format(self.angmom)
+                    f"Given convention does not support spherical contractions for the angular "
+                    f"momentum {self.angmom}"
                 )
 
             return tuple(sph_conventions[self.angmom])
@@ -227,8 +226,8 @@ def permutation_libcint(self):
 
     if molbasis.primitive_normalization != "L2":  # pragma: no cover
         raise ValueError(
-            "Only L2 normalization scheme is supported in `gbasis`. Given `IOData` instance uses "
-            "primitive normalization scheme, {}".format(molbasis.primitive_normalization)
+            f"Only L2 normalization scheme is supported in `gbasis`. Given `IOData` instance uses "
+            f"primitive normalization scheme, {molbasis.primitive_normalization}"
         )
 
     basis = []
diff --git a/pyproject.toml b/pyproject.toml
index bc2f7cc2..feff5f25 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -78,7 +78,7 @@ line-length = 100
 # RUF is Ruff-specific rules
 select = ["E", "F", "UP", "B", "I", "PGH", "PL", "RUF"]
 line-length = 100
-ignore = ["PLR2004", "PLR0913", "PLR0912", "PLW2901", "PLR0915", "RUF013"]
+ignore = ["PLR2004", "PLR0913", "PLR0912", "PLW2901", "PLR0915", "RUF013", "E741"]
 extend-exclude = ["doc/*", "doc/*/*"]
 
 [tool.pytest.ini_options]