diff --git a/.gitignore b/.gitignore
index 2ac37aea6df..8179e082ca9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,9 +1,6 @@
 __pycache__/
 .DS_Store
 pymatgen.egg-info
-dependencies/PyCifRW-3.3/PyCifRW.egg-info
-dependencies/spglib*/pyspglib.egg-info
-dependencies/spglib*/build
 *.o
 *.so
 *.pyc
@@ -26,7 +23,6 @@ setuptools*
 .cache
 .tox
 .eggs/
-gulptmp_4_1
 .coverage
 .*_cache
 # VS Code
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index bba0fa50a46..23276773f85 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -8,7 +8,7 @@ ci:
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.5.0
+    rev: v0.5.4
     hooks:
       - id: ruff
         args: [ --fix, --unsafe-fixes ]
@@ -22,7 +22,7 @@ repos:
       - id: trailing-whitespace
 
   - repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v1.10.1
+    rev: v1.11.0
     hooks:
       - id: mypy
 
@@ -65,6 +65,6 @@ repos:
         args: [ --drop-empty-cells, --keep-output ]
 
   - repo: https://github.com/RobertCraigie/pyright-python
-    rev: v1.1.369
+    rev: v1.1.373
     hooks:
       - id: pyright
diff --git a/src/pymatgen/analysis/local_env.py b/src/pymatgen/analysis/local_env.py
index d4acd5a7d2f..e5a202623d9 100644
--- a/src/pymatgen/analysis/local_env.py
+++ b/src/pymatgen/analysis/local_env.py
@@ -3048,7 +3048,7 @@ def get_order_parameters(
                                 norms[idx][j][kc] += 1
 
                         for m in range(n_neighbors):
-                            if (m != j) and (m != k) and (not flag_xaxis):
+                            if m not in {j, k} and (not flag_xaxis):
                                 tmp = max(-1.0, min(np.inner(zaxis, rij_norm[m]), 1.0))
                                 thetam = math.acos(tmp)
                                 x_two_axis_tmp = gramschmidt(rij_norm[m], zaxis)
diff --git a/src/pymatgen/analysis/magnetism/jahnteller.py b/src/pymatgen/analysis/magnetism/jahnteller.py
index 90727f19e15..ae53c08473b 100644
--- a/src/pymatgen/analysis/magnetism/jahnteller.py
+++ b/src/pymatgen/analysis/magnetism/jahnteller.py
@@ -345,7 +345,7 @@ def _get_number_of_d_electrons(species: Species) -> float:
 
         # taken from get_crystal_field_spin
         elec = species.element.full_electronic_structure
-        if len(elec) < 4 or elec[-1][1] != "s" or elec[-2][1] != "d":
+        if len(elec) < 4 or elec[-2][1] != "s" or elec[-1][1] != "d":
             raise AttributeError(f"Invalid element {species.symbol} for crystal field calculation.")
         n_electrons = int(elec[-1][2] + elec[-2][2] - species.oxi_state)  # type: ignore[operator]
         if n_electrons < 0 or n_electrons > 10:
diff --git a/src/pymatgen/core/__init__.py b/src/pymatgen/core/__init__.py
index c1f047de11c..848c251a8b4 100644
--- a/src/pymatgen/core/__init__.py
+++ b/src/pymatgen/core/__init__.py
@@ -37,9 +37,12 @@
 def _load_pmg_settings() -> dict[str, Any]:
     settings: dict[str, Any] = {}
 
+    # PMG_CONFIG_FILE takes precedence over default settings location
+    settings_file = os.getenv("PMG_CONFIG_FILE") or SETTINGS_FILE
+
     # Load .pmgrc.yaml file
     yaml = YAML()
-    for file_path in (SETTINGS_FILE, OLD_SETTINGS_FILE):
+    for file_path in (settings_file, OLD_SETTINGS_FILE):
         try:
             with open(file_path, encoding="utf-8") as yml_file:
                 settings = yaml.load(yml_file) or {}
diff --git a/src/pymatgen/core/periodic_table.py b/src/pymatgen/core/periodic_table.py
index ddc31e6bb1e..7fe90d45ce5 100644
--- a/src/pymatgen/core/periodic_table.py
+++ b/src/pymatgen/core/periodic_table.py
@@ -33,6 +33,28 @@
 
 _pt_row_sizes = (2, 8, 8, 18, 18, 32, 32)
 
+_madelung = [
+    (1, "s"),
+    (2, "s"),
+    (2, "p"),
+    (3, "s"),
+    (3, "p"),
+    (4, "s"),
+    (3, "d"),
+    (4, "p"),
+    (5, "s"),
+    (4, "d"),
+    (5, "p"),
+    (6, "s"),
+    (4, "f"),
+    (5, "d"),
+    (6, "p"),
+    (7, "s"),
+    (5, "f"),
+    (6, "d"),
+    (7, "p"),
+]
+
 
 @functools.total_ordering
 @unique
@@ -422,11 +444,12 @@ def icsd_oxidation_states(self) -> tuple[int, ...]:
     @property
     def full_electronic_structure(self) -> list[tuple[int, str, int]]:
         """Full electronic structure as list of tuples, in order of increasing
-        principal (n) and angular momentum (l)  quantum numbers.
+        energy level (according to the Madelung rule). Therefore, the final
+        element in the list gives the electronic structure of the valence shell.
 
         For example, the electronic structure for Fe is represented as:
         [(1, "s", 2), (2, "s", 2), (2, "p", 6), (3, "s", 2), (3, "p", 6),
-        (3, "d", 6), (4, "s", 2)].
+        (4, "s", 2), (3, "d", 6)].
 
         References:
             Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team (2023). NIST
@@ -445,7 +468,13 @@ def parse_orbital(orb_str):
         if data[0][0] == "[":
             sym = data[0].replace("[", "").replace("]", "")
             data = list(Element(sym).full_electronic_structure) + data[1:]
-        return data
+        # sort the final electronic structure by increasing energy level
+        return sorted(data, key=lambda x: _madelung.index((x[0], x[1])))
+
+    @property
+    def n_electrons(self) -> int:
+        """Total number of electrons in the Element."""
+        return sum([t[-1] for t in self.full_electronic_structure])
 
     @property
     def valence(self) -> tuple[int | np.nan, int]:
@@ -1117,7 +1146,8 @@ def electronic_structure(self) -> str:
     @property
     def full_electronic_structure(self) -> list[tuple[int, str, int]]:
         """Full electronic structure as list of tuples, in order of increasing
-        principal (n) and angular momentum (l)  quantum numbers.
+        energy level (according to the Madelung rule). Therefore, the final
+        element in the list gives the electronic structure of the valence shell.
 
         For example, the electronic structure for Fe+2 is represented as:
         [(1, "s", 2), (2, "s", 2), (2, "p", 6), (3, "s", 2), (3, "p", 6),
@@ -1140,7 +1170,15 @@ def parse_orbital(orb_str):
         if data[0][0] == "[":
             sym = data[0].replace("[", "").replace("]", "")
             data = list(Element(sym).full_electronic_structure) + data[1:]
-        return data
+        # sort the final electronic structure by increasing energy level
+        return sorted(data, key=lambda x: _madelung.index((x[0], x[1])))
+
+    # NOTE - copied exactly from Element. Refactoring / inheritance may improve
+    # robustness
+    @property
+    def n_electrons(self) -> int:
+        """Total number of electrons in the Species."""
+        return sum([t[-1] for t in self.full_electronic_structure])
 
     # NOTE - copied exactly from Element. Refactoring / inheritance may improve
     # robustness
@@ -1319,7 +1357,7 @@ def get_crystal_field_spin(
             raise ValueError("Invalid coordination or spin config")
 
         elec = self.element.full_electronic_structure
-        if len(elec) < 4 or elec[-1][1] != "s" or elec[-2][1] != "d":
+        if len(elec) < 4 or elec[-2][1] != "s" or elec[-1][1] != "d":
             raise AttributeError(f"Invalid element {self.symbol} for crystal field calculation")
 
         assert self.oxi_state is not None
diff --git a/src/pymatgen/electronic_structure/boltztrap.py b/src/pymatgen/electronic_structure/boltztrap.py
index 980ccb4f1a7..7f9e60f7db7 100644
--- a/src/pymatgen/electronic_structure/boltztrap.py
+++ b/src/pymatgen/electronic_structure/boltztrap.py
@@ -1,16 +1,13 @@
-"""This module provides classes to run and analyze boltztrap on pymatgen band
-structure objects. Boltztrap is a software interpolating band structures and
-computing materials properties from this band structure using Boltzmann
-semi-classical transport theory.
+"""This module provides classes to run and analyze BoltzTraP on pymatgen band
+structure objects. BoltzTraP is a software developed by Georg Madsen to
+interpolate band structures and compute materials properties from this
+band structure using Boltzmann semi-classical transport theory.
 
-Boltztrap has been developed by Georg Madsen.
-
-http://www.icams.de/content/research/software-development/boltztrap/
+https://www.tuwien.at/en/tch/tc/theoretical-materials-chemistry/boltztrap
 
 You need version 1.2.3 or higher
 
-References are:
-
+References:
     Madsen, G. K. H., and Singh, D. J. (2006).
     BoltzTraP. A code for calculating band-structure dependent quantities.
     Computer Physics Communications, 175, 67-71
@@ -60,13 +57,13 @@
 
 
 class BoltztrapRunner(MSONable):
-    """This class is used to run Boltztrap on a band structure object."""
+    """This class is used to run BoltzTraP on a band structure object."""
 
     @requires(
         which("x_trans"),
         "BoltztrapRunner requires the executables 'x_trans' to be in PATH. Please download "
-        "Boltztrap at http://www.icams.de/content/research/software-development/boltztrap/ "
-        "and follow the instructions in the README to compile Bolztrap accordingly. "
+        "BoltzTraP at https://www.tuwien.at/en/tch/tc/theoretical-materials-chemistry/boltztrap "
+        "and follow the instructions in the README to compile BoltzTraP accordingly. "
         "Then add x_trans to your path",
     )
     def __init__(
@@ -144,7 +141,7 @@ def __init__(
                 electron occupations. If the band structure comes from a soc
                 computation, you should set soc to True (default False)
             doping:
-                the fixed doping levels you want to compute. Boltztrap provides
+                the fixed doping levels you want to compute. BoltzTraP provides
                 both transport values depending on electron chemical potential
                 (fermi energy) and for a series of fixed carrier
                 concentrations. By default, this is set to 1e16 to 1e22 in
@@ -734,7 +731,7 @@ def __init__(
         bz_kpoints=None,
         fermi_surface_data=None,
     ) -> None:
-        """Constructor taking directly all the data generated by Boltztrap. You
+        """Constructor taking directly all the data generated by BoltzTraP. You
         won't probably use it directly but instead use the from_files and
         from_dict methods.
 
@@ -760,7 +757,7 @@ def __init__(
                 each Fermi level in mu_steps]}
                 The units are m^3/C
             doping: The different doping levels that have been given to
-                Boltztrap. The format is {'p':[],'n':[]} with an array of
+                BoltzTraP. The format is {'p':[],'n':[]} with an array of
                 doping levels. The units are cm^-3
             mu_doping: Gives the electron chemical potential (or Fermi level)
                 for a given set of doping.
@@ -802,7 +799,7 @@ def __init__(
             intrans: a dictionary of inputs e.g. {"scissor": 0.0}
             carrier_conc: The concentration of carriers in electron (or hole)
                 per unit cell
-            dos: The dos computed by Boltztrap given as a pymatgen Dos object
+            dos: The dos computed by BoltzTraP given as a pymatgen Dos object
             dos_partial: Data for the partial DOS projected on sites and
                 orbitals
             vol: Volume of the unit cell in angstrom cube (A^3)
@@ -837,13 +834,13 @@ def __init__(
         self.fermi_surface_data = fermi_surface_data
 
     def get_symm_bands(self, structure: Structure, efermi, kpt_line=None, labels_dict=None):
-        """Useful to read bands from Boltztrap output and get a BandStructureSymmLine object
+        """Useful to read bands from BoltzTraP output and get a BandStructureSymmLine object
         comparable with that one from a DFT calculation (if the same kpt_line is
         provided). Default kpt_line and labels_dict is the standard path of high symmetry
         k-point for the specified structure. They could be extracted from the
         BandStructureSymmLine object that you want to compare with. efermi variable must
         be specified to create the BandStructureSymmLine object (usually it comes from DFT
-        or Boltztrap calc).
+        or BoltzTraP calc).
         """
         try:
             if kpt_line is None:
@@ -1640,7 +1637,7 @@ def get_carrier_concentration(self):
 
     def get_hall_carrier_concentration(self):
         """Get the Hall carrier concentration (in cm^-3). This is the trace of
-        the Hall tensor (see Boltztrap source code) Hall carrier concentration
+        the Hall tensor (see BoltzTraP source code) Hall carrier concentration
         are not always exactly the same than carrier concentration.
 
         Returns:
@@ -1770,7 +1767,7 @@ def parse_intrans(path_dir):
             path_dir: (str) dir containing the boltztrap.intrans file
 
         Returns:
-            dict: various inputs that had been used in the Boltztrap run.
+            dict: various inputs that had been used in the BoltzTraP run.
         """
         intrans = {}
         with open(f"{path_dir}/boltztrap.intrans") as file:
diff --git a/src/pymatgen/io/aims/parsers.py b/src/pymatgen/io/aims/parsers.py
index 3021c49279c..0a70723e57e 100644
--- a/src/pymatgen/io/aims/parsers.py
+++ b/src/pymatgen/io/aims/parsers.py
@@ -330,7 +330,7 @@ def _parse_k_points(self) -> None:
 
         line_start = self.reverse_search_for(["| K-points in task"])
         line_end = self.reverse_search_for(["| k-point:"])
-        if (line_start == LINE_NOT_FOUND) or (line_end == LINE_NOT_FOUND) or (line_end - line_start != n_kpts):
+        if LINE_NOT_FOUND in {line_start, line_end} or (line_end - line_start != n_kpts):
             self._cache.update(
                 {
                     "k_points": None,
diff --git a/src/pymatgen/io/vasp/inputs.py b/src/pymatgen/io/vasp/inputs.py
index 25572a2b1ab..655396e8e0c 100644
--- a/src/pymatgen/io/vasp/inputs.py
+++ b/src/pymatgen/io/vasp/inputs.py
@@ -2856,7 +2856,7 @@ def from_directory(
                 dict of {filename: Object type}. Objects must have
                 from_file method.
         """
-        sub_dct = {}
+        sub_dct: dict[str, Any] = {}
         for fname, ftype in (
             ("INCAR", Incar),
             ("KPOINTS", Kpoints),
diff --git a/src/pymatgen/io/vasp/outputs.py b/src/pymatgen/io/vasp/outputs.py
index 9df9d6dadc9..9466ca1b7d2 100644
--- a/src/pymatgen/io/vasp/outputs.py
+++ b/src/pymatgen/io/vasp/outputs.py
@@ -1703,10 +1703,8 @@ def __init__(
                 tag = elem.tag
                 if event == "start":
                     # The start event tells us when we have entered blocks
-                    if (
-                        tag == "eigenvalues_kpoints_opt"
-                        or tag == "projected_kpoints_opt"
-                        or (tag == "dos" and elem.attrib.get("comment") == "kpoints_opt")
+                    if tag in {"eigenvalues_kpoints_opt", "projected_kpoints_opt"} or (
+                        tag == "dos" and elem.attrib.get("comment") == "kpoints_opt"
                     ):
                         in_kpoints_opt = True
                 elif not parsed_header:
@@ -3891,31 +3889,31 @@ def __init__(self, filename: PathLike) -> None:
             self.data = data
             self.phase_factors = phase_factors
 
-    def get_projection_on_elements(self, structure: Structure) -> dict[Spin, list]:
+    def get_projection_on_elements(self, structure: Structure) -> dict[Spin, list[list[dict[str, float]]]]:
         """Get a dict of projections on elements.
 
         Args:
             structure (Structure): Input structure.
 
         Returns:
-            A dict as {Spin.up: [k index][b index][{Element: values}]].
+            A dict as {Spin: [band index][kpoint index][{Element: values}]].
         """
         assert self.data is not None, "Data cannot be None."
         assert self.nkpoints is not None
         assert self.nbands is not None
         assert self.nions is not None
 
-        dico: dict[Spin, list] = {}
+        elem_proj: dict[Spin, list] = {}
         for spin in self.data:
-            dico[spin] = [[defaultdict(float) for _ in range(self.nkpoints)] for _ in range(self.nbands)]
+            elem_proj[spin] = [[defaultdict(float) for _ in range(self.nkpoints)] for _ in range(self.nbands)]
 
         for iat in range(self.nions):
             name = structure.species[iat].symbol
             for spin, data in self.data.items():
                 for kpoint, band in itertools.product(range(self.nkpoints), range(self.nbands)):
-                    dico[spin][band][kpoint][name] += np.sum(data[kpoint, band, iat, :])
+                    elem_proj[spin][band][kpoint][name] += np.sum(data[kpoint, band, iat, :])
 
-        return dico
+        return elem_proj
 
     def get_occupation(self, atom_index: int, orbital: str) -> dict:
         """Get the occupation for a particular orbital of a particular atom.
diff --git a/src/pymatgen/transformations/advanced_transformations.py b/src/pymatgen/transformations/advanced_transformations.py
index 57b89c43efe..c6107edc562 100644
--- a/src/pymatgen/transformations/advanced_transformations.py
+++ b/src/pymatgen/transformations/advanced_transformations.py
@@ -1437,29 +1437,42 @@ def __init__(
         min_atoms: int | None = None,
         max_atoms: int | None = None,
         min_length: float = 15.0,
+        max_length: float | None = None,
         force_diagonal: bool = False,
         force_90_degrees: bool = False,
+        allow_orthorhombic: bool = False,
         angle_tolerance: float = 1e-3,
+        step_size: float = 0.1,
     ):
         """
         Args:
             max_atoms: Maximum number of atoms allowed in the supercell.
             min_atoms: Minimum number of atoms allowed in the supercell.
             min_length: Minimum length of the smallest supercell lattice vector.
+            max_length: Maximum length of the larger supercell lattice vector.
             force_diagonal: If True, return a transformation with a diagonal
                 transformation matrix.
             force_90_degrees: If True, return a transformation for a supercell
                 with 90 degree angles (if possible). To avoid long run times,
-                please use max_atoms
+                please use max_atoms or max_length
+            allow_orthorhombic: Instead of a cubic cell, also orthorhombic cells
+                are allowed. max_length is required for this option.
             angle_tolerance: tolerance to determine the 90 degree angles.
+            step_size (float): step_size which is used to increase the supercell.
+                If allow_orthorhombic and force_90_degrees is both set to True,
+                the chosen step_size will be automatically multiplied by 5 to
+                prevent a too long search for the possible supercell.
         """
         self.min_atoms = min_atoms or -np.inf
         self.max_atoms = max_atoms or np.inf
         self.min_length = min_length
+        self.max_length = max_length
         self.force_diagonal = force_diagonal
         self.force_90_degrees = force_90_degrees
+        self.allow_orthorhombic = allow_orthorhombic
         self.angle_tolerance = angle_tolerance
         self.transformation_matrix = None
+        self.step_size = step_size
 
     def apply_transformation(self, structure: Structure) -> Structure:
         """The algorithm solves for a transformation matrix that makes the
@@ -1478,8 +1491,8 @@ def apply_transformation(self, structure: Structure) -> Structure:
         """
         lat_vecs = structure.lattice.matrix
 
-        # boolean for if a sufficiently large supercell has been created
-        sc_not_found = True
+        if self.max_length is None and self.allow_orthorhombic:
+            raise AttributeError("max_length is required for orthorhombic cells")
 
         if self.force_diagonal:
             scale = self.min_length / np.array(structure.lattice.abc)
@@ -1487,65 +1500,120 @@ def apply_transformation(self, structure: Structure) -> Structure:
             st = SupercellTransformation(self.transformation_matrix)
             return st.apply_transformation(structure)
 
-        # target_threshold is used as the desired cubic side lengths
-        target_sc_size = self.min_length
-        while sc_not_found:
-            target_sc_lat_vecs = np.eye(3, 3) * target_sc_size
-            self.transformation_matrix = target_sc_lat_vecs @ np.linalg.inv(lat_vecs)
+        if not self.allow_orthorhombic:
+            # boolean for if a sufficiently large supercell has been created
+            sc_not_found = True
 
-            # round the entries of T and force T to be non-singular
-            self.transformation_matrix = _round_and_make_arr_singular(  # type: ignore[assignment]
-                self.transformation_matrix  # type: ignore[arg-type]
+            # target_threshold is used as the desired cubic side lengths
+            target_sc_size = self.min_length
+            while sc_not_found:
+                target_sc_lat_vecs = np.eye(3, 3) * target_sc_size
+                length_vecs, n_atoms, superstructure, self.transformation_matrix = self.get_possible_supercell(
+                    lat_vecs, structure, target_sc_lat_vecs
+                )
+                # Check if constraints are satisfied
+                if self.check_constraints(length_vecs=length_vecs, n_atoms=n_atoms, superstructure=superstructure):
+                    return superstructure
+
+                # Increase threshold until proposed supercell meets requirements
+                target_sc_size += self.step_size
+                self.check_exceptions(length_vecs, n_atoms)
+
+            raise AttributeError("Unable to find cubic supercell")
+
+        if self.force_90_degrees:
+            # prevent a too long search for the supercell
+            self.step_size *= 5
+
+        combined_list = [
+            [size_a, size_b, size_c]
+            for size_a in np.arange(self.min_length, self.max_length, self.step_size)
+            for size_b in np.arange(self.min_length, self.max_length, self.step_size)
+            for size_c in np.arange(self.min_length, self.max_length, self.step_size)
+        ]
+        combined_list = sorted(combined_list, key=sum)
+
+        for size_a, size_b, size_c in combined_list:
+            target_sc_lat_vecs = np.array([[size_a, 0, 0], [0, size_b, 0], [0, 0, size_c]])
+            length_vecs, n_atoms, superstructure, self.transformation_matrix = self.get_possible_supercell(
+                lat_vecs, structure, target_sc_lat_vecs
             )
+            # Check if constraints are satisfied
+            if self.check_constraints(length_vecs=length_vecs, n_atoms=n_atoms, superstructure=superstructure):
+                return superstructure
+
+            self.check_exceptions(length_vecs, n_atoms)
+        raise AttributeError("Unable to find orthorhombic supercell")
 
-            proposed_sc_lat_vecs = self.transformation_matrix @ lat_vecs
-
-            # Find the shortest dimension length and direction
-            a = proposed_sc_lat_vecs[0]
-            b = proposed_sc_lat_vecs[1]
-            c = proposed_sc_lat_vecs[2]
-
-            length1_vec = c - _proj(c, a)  # a-c plane
-            length2_vec = a - _proj(a, c)
-            length3_vec = b - _proj(b, a)  # b-a plane
-            length4_vec = a - _proj(a, b)
-            length5_vec = b - _proj(b, c)  # b-c plane
-            length6_vec = c - _proj(c, b)
-            length_vecs = np.array(
-                [
-                    length1_vec,
-                    length2_vec,
-                    length3_vec,
-                    length4_vec,
-                    length5_vec,
-                    length6_vec,
-                ]
+    def check_exceptions(self, length_vecs, n_atoms):
+        """Check supercell exceptions."""
+        if n_atoms > self.max_atoms:
+            raise AttributeError(
+                "While trying to solve for the supercell, the max "
+                "number of atoms was exceeded. Try lowering the number"
+                "of nearest neighbor distances."
             )
+        if self.max_length is not None and np.max(np.linalg.norm(length_vecs, axis=1)) >= self.max_length:
+            raise AttributeError("While trying to solve for the supercell, the max length was exceeded.")
 
-            # Get number of atoms
-            st = SupercellTransformation(self.transformation_matrix)
-            superstructure = st.apply_transformation(structure)
-            n_atoms = len(superstructure)
+    def check_constraints(self, length_vecs, n_atoms, superstructure):
+        """
+        Check if the supercell constraints are met.
 
-            # Check if constraints are satisfied
-            if (
+        Returns:
+            bool
+
+        """
+        return bool(
+            (
                 np.min(np.linalg.norm(length_vecs, axis=1)) >= self.min_length
                 and self.min_atoms <= n_atoms <= self.max_atoms
-            ) and (
+            )
+            and (
                 not self.force_90_degrees
                 or np.all(np.absolute(np.array(superstructure.lattice.angles) - 90) < self.angle_tolerance)
-            ):
-                return superstructure
+            )
+        )
 
-            # Increase threshold until proposed supercell meets requirements
-            target_sc_size += 0.1
-            if n_atoms > self.max_atoms:
-                raise AttributeError(
-                    "While trying to solve for the supercell, the max "
-                    "number of atoms was exceeded. Try lowering the number"
-                    "of nearest neighbor distances."
-                )
-        raise AttributeError("Unable to find cubic supercell")
+    @staticmethod
+    def get_possible_supercell(lat_vecs, structure, target_sc_lat_vecs):
+        """
+        Get the supercell possible with the set conditions.
+
+        Returns:
+            length_vecs, n_atoms, superstructure, transformation_matrix
+        """
+        transformation_matrix = target_sc_lat_vecs @ np.linalg.inv(lat_vecs)
+        # round the entries of T and force T to be non-singular
+        transformation_matrix = _round_and_make_arr_singular(  # type: ignore[assignment]
+            transformation_matrix  # type: ignore[arg-type]
+        )
+        proposed_sc_lat_vecs = transformation_matrix @ lat_vecs
+        # Find the shortest dimension length and direction
+        a = proposed_sc_lat_vecs[0]
+        b = proposed_sc_lat_vecs[1]
+        c = proposed_sc_lat_vecs[2]
+        length1_vec = c - _proj(c, a)  # a-c plane
+        length2_vec = a - _proj(a, c)
+        length3_vec = b - _proj(b, a)  # b-a plane
+        length4_vec = a - _proj(a, b)
+        length5_vec = b - _proj(b, c)  # b-c plane
+        length6_vec = c - _proj(c, b)
+        length_vecs = np.array(
+            [
+                length1_vec,
+                length2_vec,
+                length3_vec,
+                length4_vec,
+                length5_vec,
+                length6_vec,
+            ]
+        )
+        # Get number of atoms
+        st = SupercellTransformation(transformation_matrix)
+        superstructure = st.apply_transformation(structure)
+        n_atoms = len(superstructure)
+        return length_vecs, n_atoms, superstructure, transformation_matrix
 
 
 class AddAdsorbateTransformation(AbstractTransformation):
diff --git a/src/pymatgen/util/testing/__init__.py b/src/pymatgen/util/testing/__init__.py
index 211690effa9..04aea7094fa 100644
--- a/src/pymatgen/util/testing/__init__.py
+++ b/src/pymatgen/util/testing/__init__.py
@@ -49,7 +49,7 @@ def _tmp_dir(self, tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None:
     @classmethod
     def get_structure(cls, name: str) -> Structure:
         """
-        Lazily load a structure from pymatgen/util/testing/structures.
+        Lazily load a structure from pymatgen/util/structures.
 
         Args:
             name (str): Name of structure file.
diff --git a/tests/core/test_periodic_table.py b/tests/core/test_periodic_table.py
index f92335d5266..586a6a513bf 100644
--- a/tests/core/test_periodic_table.py
+++ b/tests/core/test_periodic_table.py
@@ -74,8 +74,8 @@ def test_full_electronic_structure(self):
                 (2, "p", 6),
                 (3, "s", 2),
                 (3, "p", 6),
-                (3, "d", 6),
                 (4, "s", 2),
+                (3, "d", 6),
             ],
             "Li": [(1, "s", 2), (2, "s", 1)],
             "U": [
@@ -84,19 +84,19 @@ def test_full_electronic_structure(self):
                 (2, "p", 6),
                 (3, "s", 2),
                 (3, "p", 6),
-                (3, "d", 10),
                 (4, "s", 2),
+                (3, "d", 10),
                 (4, "p", 6),
-                (4, "d", 10),
                 (5, "s", 2),
+                (4, "d", 10),
                 (5, "p", 6),
+                (6, "s", 2),
                 (4, "f", 14),
                 (5, "d", 10),
-                (6, "s", 2),
                 (6, "p", 6),
+                (7, "s", 2),
                 (5, "f", 3),
                 (6, "d", 1),
-                (7, "s", 2),
             ],
         }
         for k, v in cases.items():
@@ -169,6 +169,11 @@ def test_from_row_and_group(self):
         for k, v in cases.items():
             assert ElementBase.from_row_and_group(v[0], v[1]) == Element(k)
 
+    def test_n_electrons(self):
+        cases = {"O": 8, "Fe": 26, "Li": 3, "Be": 4}
+        for k, v in cases.items():
+            assert Element(k).n_electrons == v
+
     def test_valence(self):
         cases = {"O": (1, 4), "Fe": (2, 6), "Li": (0, 1), "Be": (0, 2)}
         for k, v in cases.items():
@@ -602,18 +607,20 @@ def test_sort(self):
 
     def test_species_electronic_structure(self):
         assert Species("Fe", 0).electronic_structure == "[Ar].3d6.4s2"
+        assert Species("Fe", 0).n_electrons == 26
         assert Species("Fe", 0).full_electronic_structure == [
             (1, "s", 2),
             (2, "s", 2),
             (2, "p", 6),
             (3, "s", 2),
             (3, "p", 6),
-            (3, "d", 6),
             (4, "s", 2),
+            (3, "d", 6),
         ]
         assert Species("Fe", 0).valence == (2, 6)
 
         assert Species("Fe", 2).electronic_structure == "[Ar].3d6"
+        assert Species("Fe", 2).n_electrons == 24
         assert Species("Fe", 2).full_electronic_structure == [
             (1, "s", 2),
             (2, "s", 2),
@@ -625,6 +632,7 @@ def test_species_electronic_structure(self):
         assert Species("Fe", 2).valence == (2, 6)
 
         assert Species("Fe", 3).electronic_structure == "[Ar].3d5"
+        assert Species("Fe", 3).n_electrons == 23
         assert Species("Fe", 3).full_electronic_structure == [
             (1, "s", 2),
             (2, "s", 2),
@@ -635,12 +643,36 @@ def test_species_electronic_structure(self):
         ]
         assert Species("Fe", 3).valence == (2, 5)
 
+        assert Species("Th", 4).electronic_structure == "[Hg].6p6"
+        assert Species("Th", 4).full_electronic_structure == [
+            (1, "s", 2),
+            (2, "s", 2),
+            (2, "p", 6),
+            (3, "s", 2),
+            (3, "p", 6),
+            (4, "s", 2),
+            (3, "d", 10),
+            (4, "p", 6),
+            (5, "s", 2),
+            (4, "d", 10),
+            (5, "p", 6),
+            (6, "s", 2),
+            (4, "f", 14),
+            (5, "d", 10),
+            (6, "p", 6),
+        ]
+        assert Species("Th", 4).valence == (1, 6)
+
         assert Species("Li", 1).electronic_structure == "1s2"
+        assert Species("Li", 1).n_electrons == 2
         # alkali metals, all p
         for el in ["Na", "K", "Rb", "Cs"]:
             assert Species(el, 1).electronic_structure.split(".")[-1][1::] == "p6", f"Failure for {el} +1"
         for el in ["Ca", "Mg", "Ba", "Sr"]:
             assert Species(el, 2).electronic_structure.split(".")[-1][1::] == "p6", f"Failure for {el} +2"
+        # valence shell should be f (l=3) for all lanthanide ions except La+3 and Lu+3
+        for el in ["Ce", "Nd", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu"]:
+            assert Species(el, 3).valence[0] == 3, f"Failure for {el} +3"
 
         for el in Element:
             for ox in el.common_oxidation_states:
diff --git a/tests/core/test_settings.py b/tests/core/test_settings.py
index 39107becc90..6311909f6b1 100644
--- a/tests/core/test_settings.py
+++ b/tests/core/test_settings.py
@@ -48,3 +48,34 @@ def test_load_settings(tmp_path: Path, monkeypatch: MonkeyPatch) -> None:
     # should return empty dict if file is invalid
     settings_file.write_text("---")
     assert _load_pmg_settings() == {}
+
+
+def test_env_var_pmg_config_file(tmp_path: Path, monkeypatch: MonkeyPatch) -> None:
+    custom_config_file = tmp_path / "custom_config.yaml"
+    custom_config_file.write_text("PMG_CUSTOM_SETTING: custom_value")
+
+    with monkeypatch.context() as ctx:
+        ctx.setenv("PMG_CONFIG_FILE", str(custom_config_file))
+        settings = _load_pmg_settings()
+        assert "PMG_CUSTOM_SETTING" in settings
+        assert settings["PMG_CUSTOM_SETTING"] == "custom_value"
+
+    # Test that PMG_CONFIG_FILE takes precedence over the default location
+    settings_file = tmp_path / ".pmgrc.yaml"
+    monkeypatch.setattr("pymatgen.core.SETTINGS_FILE", settings_file)
+    settings_file.write_text("PMG_DEFAULT_SETTING: default_value")
+    custom_config_file.write_text("PMG_CUSTOM_SETTING: custom_value")
+
+    with monkeypatch.context() as ctx:
+        ctx.setenv("PMG_CONFIG_FILE", str(custom_config_file))
+        settings = _load_pmg_settings()
+        assert "PMG_CUSTOM_SETTING" in settings
+        assert "PMG_DEFAULT_SETTING" not in settings
+        assert settings["PMG_CUSTOM_SETTING"] == "custom_value"
+
+    # Test that env vars still take precedence over the values specified in PMG_CONFIG_FILE
+    with monkeypatch.context() as ctx:
+        ctx.setenv("PMG_CONFIG_FILE", str(custom_config_file))
+        ctx.setenv("PMG_CUSTOM_SETTING", "env_value")
+        settings = _load_pmg_settings()
+        assert settings["PMG_CUSTOM_SETTING"] == "env_value"
diff --git a/tests/io/vasp/test_inputs.py b/tests/io/vasp/test_inputs.py
index aba859bea1f..4ce0c20e95f 100644
--- a/tests/io/vasp/test_inputs.py
+++ b/tests/io/vasp/test_inputs.py
@@ -1099,8 +1099,8 @@ def test_nelectrons(self):
         assert self.psingle_Fe.nelectrons == 8
 
     def test_electron_config(self):
-        assert self.psingle_Mn_pv.electron_configuration == [(4, "s", 2), (3, "d", 5), (3, "p", 6)]
-        assert self.psingle_Fe.electron_configuration == [(4, "s", 2), (3, "d", 6)]
+        assert self.psingle_Mn_pv.electron_configuration == [(3, "d", 5), (4, "s", 2), (3, "p", 6)]
+        assert self.psingle_Fe.electron_configuration == [(3, "d", 6), (4, "s", 2)]
 
     def test_attributes(self):
         for key, val in self.Mn_pv_attrs.items():
diff --git a/tests/transformations/test_advanced_transformations.py b/tests/transformations/test_advanced_transformations.py
index ea0feed09b8..653c098bb19 100644
--- a/tests/transformations/test_advanced_transformations.py
+++ b/tests/transformations/test_advanced_transformations.py
@@ -692,7 +692,7 @@ def test_monte_carlo(self):
 
 
 class TestCubicSupercellTransformation(PymatgenTest):
-    def test_apply_transformation(self):
+    def test_apply_transformation_cubic_supercell(self):
         structure = self.get_structure("TlBiSe2")
         min_atoms = 100
         max_atoms = 1000
@@ -756,6 +756,84 @@ def test_apply_transformation(self):
         transformed_structure = supercell_generator.apply_transformation(structure)
         assert_allclose(list(transformed_structure.lattice.angles), [90.0, 90.0, 90.0])
 
+    def test_apply_transformation_orthorhombic_supercell(self):
+        structure = self.get_structure("Li3V2(PO4)3")
+        min_atoms = 100
+        max_atoms = 400
+
+        supercell_generator_cubic = CubicSupercellTransformation(
+            min_atoms=min_atoms,
+            max_atoms=max_atoms,
+            min_length=10.0,
+            force_90_degrees=False,
+            allow_orthorhombic=False,
+            max_length=25,
+        )
+
+        transformed_cubic = supercell_generator_cubic.apply_transformation(structure)
+
+        supercell_generator_orthorhombic = CubicSupercellTransformation(
+            min_atoms=min_atoms,
+            max_atoms=max_atoms,
+            min_length=10.0,
+            force_90_degrees=False,
+            allow_orthorhombic=True,
+            max_length=25,
+        )
+
+        transformed_orthorhombic = supercell_generator_orthorhombic.apply_transformation(structure)
+
+        assert_array_equal(
+            supercell_generator_orthorhombic.transformation_matrix,
+            np.array([[0, -2, 1], [-2, 0, 0], [0, 0, -2]]),
+        )
+
+        # make sure that the orthorhombic supercell is different from the cubic cell
+        assert not np.array_equal(
+            supercell_generator_cubic.transformation_matrix, supercell_generator_orthorhombic.transformation_matrix
+        )
+        assert transformed_cubic.lattice.angles != transformed_orthorhombic.lattice.angles
+        assert transformed_orthorhombic.lattice.abc != transformed_cubic.lattice.abc
+
+        structure = self.get_structure("Si")
+        min_atoms = 100
+        max_atoms = 400
+
+        supercell_generator_cubic = CubicSupercellTransformation(
+            min_atoms=min_atoms,
+            max_atoms=max_atoms,
+            min_length=10.0,
+            force_90_degrees=True,
+            allow_orthorhombic=False,
+            max_length=25,
+        )
+
+        transformed_cubic = supercell_generator_cubic.apply_transformation(structure)
+
+        supercell_generator_orthorhombic = CubicSupercellTransformation(
+            min_atoms=min_atoms,
+            max_atoms=max_atoms,
+            min_length=10.0,
+            force_90_degrees=True,
+            allow_orthorhombic=True,
+            max_length=25,
+        )
+
+        transformed_orthorhombic = supercell_generator_orthorhombic.apply_transformation(structure)
+
+        assert_array_equal(
+            supercell_generator_orthorhombic.transformation_matrix,
+            np.array([[3, 0, 0], [-2, 4, 0], [-2, 4, 6]]),
+        )
+
+        # make sure that the orthorhombic supercell is different from the cubic cell
+        assert not np.array_equal(
+            supercell_generator_cubic.transformation_matrix, supercell_generator_orthorhombic.transformation_matrix
+        )
+        assert transformed_orthorhombic.lattice.abc != transformed_cubic.lattice.abc
+        # only angels are expected to be the same because of force_90_degrees = True
+        assert transformed_cubic.lattice.angles == transformed_orthorhombic.lattice.angles
+
 
 class TestAddAdsorbateTransformation(PymatgenTest):
     def test_apply_transformation(self):