diff --git a/src/koopmans/calculators/_projwfc.py b/src/koopmans/calculators/_projwfc.py
index 1338cdbd1..81181bf07 100644
--- a/src/koopmans/calculators/_projwfc.py
+++ b/src/koopmans/calculators/_projwfc.py
@@ -75,7 +75,8 @@ def generate_dos(self) -> GridDOSCollection:
             # The filename does not encode the principal quantum number n. In order to recover this number, we compare
             # the reported angular momentum quantum number l against the list of expected orbitals, and infer n
             # assuming only that the file corresponding to nl will come before (n+1)l
-            orbitals = copy.copy(self._expected_orbitals[atom.symbol])
+            label = atom.symbol + str(atom.tag) if atom.tag > 0 else atom.symbol
+            orbitals = copy.copy(self._expected_orbitals[label])
             for filename in filenames:
                 # Infer l from the filename
                 subshell = filename.name[-2]
diff --git a/src/koopmans/pseudopotentials.py b/src/koopmans/pseudopotentials.py
index 45b1c2e9d..b56046190 100644
--- a/src/koopmans/pseudopotentials.py
+++ b/src/koopmans/pseudopotentials.py
@@ -160,14 +160,15 @@ def expected_subshells(atoms: Atoms, pseudopotentials: Dict[str, str],
 
     expected_orbitals = {}
     for atom in atoms:
-        if atom.symbol in expected_orbitals:
+        label = atom.symbol + str(atom.tag) if atom.tag > 0 else atom.symbol
+        if label in expected_orbitals:
             continue
-        pseudo_file = pseudopotentials[atom.symbol]
+        pseudo_file = pseudopotentials[label]
         z_core = atom.number - valence_from_pseudo(pseudo_file, pseudo_dir)
         if z_core in z_core_to_first_orbital:
             first_orbital = z_core_to_first_orbital[z_core]
         else:
             raise ValueError(f'Failed to identify the subshells of the valence of {pseudo_file}')
         all_orbitals = list(z_core_to_first_orbital.values()) + ['5d', '6p', '6d']
-        expected_orbitals[atom.symbol] = sorted(all_orbitals[all_orbitals.index(first_orbital):])
+        expected_orbitals[label] = sorted(all_orbitals[all_orbitals.index(first_orbital):])
     return expected_orbitals