diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py
index 85c45d18..e768caa8 100644
--- a/src/tequila/quantumchemistry/chemistry_tools.py
+++ b/src/tequila/quantumchemistry/chemistry_tools.py
@@ -302,13 +302,13 @@ def get_number_of_core_electrons(self):
             if n > 10:
                 result += 10 - 2
             if n > 18:
-                result += 18 - 10 - 2
+                result += 18 - 10
             if n > 36:
-                result += 36 - 18 - 10 - 2
+                result += 36 - 18
             if n > 54:
-                result += 54 - 36 - 18 - 10 - 2
+                result += 54 - 36
             if n > 86:
-                result += 86 - 54 - 36 - 18 - 10 - 2
+                result += 86 - 54
         return result
 
     @property
diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 6151572c..5466e217 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -685,7 +685,7 @@ def orthonormalize_basis_orbitals(self):
         # backward compatibility
         return self.use_native_orbitals()
 
-    def use_native_orbitals(self, inplace=False, core: list = [], *args, **kwargs):
+    def use_native_orbitals(self, inplace=False, core: list = None, *args, **kwargs):
         """
         Parameters
         ----------
@@ -710,118 +710,116 @@ def use_native_orbitals(self, inplace=False, core: list = [], *args, **kwargs):
         New molecule in the native (orthonormalized) basis given
         e.g. for standard basis sets the orbitals are orthonormalized Gaussian Basis Functions
         """
-        c = copy.deepcopy(self.integral_manager.orbital_coefficients)
-        s = self.integral_manager.overlap_integrals
-        d = self.integral_manager.get_orthonormalized_orbital_coefficients()
+        c = self.integral_manager.orbital_coefficients.copy()
+        s = self.integral_manager.overlap_integrals.copy()
+        d = self.integral_manager.get_orthonormalized_orbital_coefficients().copy()
 
-        def inner(a, b, s):
-            return numpy.sum(numpy.multiply(numpy.outer(a, b), s))
+        def orthogonalize_active_space(
+            c: numpy.ndarray, s: numpy.ndarray, frozen_idx: list[int], active_idx: list[int]
+        ) -> numpy.ndarray:
+            """
+            Symmetrically orthogonalize orbital coefficients defined by colums with indices 'active_idx' while keeping untouched those defined by 'frozen_idx'.
+            c: (basis_functions, orbitals)
+            s: (basis_functions, basis_functions)
+            frozen_idx: list of orbitals to left untoched
+            active_idx: list of orbitals to orthogonalize
+            """
+            c_f = c[:, frozen_idx]
+            c_a = c[:, active_idx]
+
+            sf = c_f.T @ s @ c_f
+            cross = c_f.T @ s @ c_a
+
+            proj = c_f @ numpy.linalg.solve(sf, cross)
+            c_a_proj = c_a - proj
+
+            sa = c_a_proj.T @ s @ c_a_proj
+            e, U = numpy.linalg.eigh(sa)
+            e = numpy.maximum(e, 1e-12)
+            X = U @ numpy.diag(1.0 / numpy.sqrt(e)) @ U.T
+
+            c_new = c.copy()
+            c_new[:, active_idx] = c_a_proj @ X
+            return c_new
 
-        def orthogonalize(c, d, s):
+        def get_active(
+            c_orig: numpy.ndarray, d_orig: numpy.ndarray, s: numpy.ndarray, active_idx_c: list[int]
+        ) -> list[int]:
             """
-            :return: orthogonalized orbitals with core HF orbitals and active Orthongonalized Native orbitals.
+            Safely identifies active orbitals in d_orig by projecting them into the entire active subspace of c_orig.
+            c_orig: original orbital matrix which will be frozen (typicall HF)
+            d_orig: original orbital matrix to look for active w.r.t. c_orig (i.e. native orbital matrix or CLPO matrix before active space considerations)
+            s: overlap_integrals
+            active_idx_c: subspace from c_orig to  look for the active indices for d_orig
             """
-            ### Computing Core-Active overlap Matrix
-            # sbar_{ki} = \langle \phi_k | \varphi_i \rangle = \sum_{m,n} c_{nk}d_{mi}\langle \chi_n | \chi_m \rangle
-            # c_{nk} = HF coeffs, d_{mi} = nat orb coef s_{mn} = Atomic Overlap Matrix
-            # k \in active orbs, i \in core orbs, m,n \in basis coeffs
-            # sbar = np.einsum('nk,mi,nm->ki', c, d, s) #only works if active == to_active
-            c = c.T
-            d = d.T
-            sbar = numpy.zeros(shape=s.shape)
-            for k in active:
-                for i in core:
-                    sbar[i][to_active[k]] = inner(c[i], d[k], s)
-            ### Projecting out Core orbitals from the Native ones
-            # dbar_{ji} = d_{ji} - \sum_k sbar_{ki}c_{jk}
-            # k \in active, i \in core, j in basis coeffs
-            dbar = numpy.zeros(shape=s.shape)
-
-            for j in active:
-                dbar[to_active[j]] = d[j]
-                for i in core:
-                    temp = sbar[i][to_active[j]] * c[i]
-                    dbar[to_active[j]] -= temp
-            ### Projected-out Nat Orbs Normalization
-            for i in to_active.values():
-                norm = numpy.sqrt(numpy.sum(numpy.multiply(numpy.outer(dbar[i], dbar[i]), s.T)))
-                if not numpy.isclose(norm, 0):
-                    dbar[i] = dbar[i] / norm
-            ### Reintroducing the New Coeffs on the HF coeff matrix
-            for j in to_active.values():
-                c[j] = dbar[j]
-            ### Compute new orbital overlap matrix:
-            sprima = numpy.eye(len(c))
-            for idx, i in enumerate(to_active.values()):
-                for j in [*to_active.values()][idx:]:
-                    sprima[i][j] = inner(c[i], c[j], s)
-                    sprima[j][i] = sprima[i][j]
-            ### Symmetric orthonormalization
-            lam_s, l_s = numpy.linalg.eigh(sprima)
-            lam_s = lam_s * numpy.eye(len(lam_s))
-            lam_sqrt_inv = numpy.sqrt(numpy.linalg.inv(lam_s))
-            symm_orthog = numpy.dot(l_s, numpy.dot(lam_sqrt_inv, l_s.T))
-            return symm_orthog.dot(c).T
-
-        def get_active(core):
-            ov = numpy.zeros(shape=(len(self.integral_manager.orbitals)))
-            for i in core:
-                for j in range(len(d)):
-                    ov[j] += numpy.abs(inner(c.T[i], d.T[j], s))
-            act = []
-            for i in range(len(self.integral_manager.orbitals) - len(core)):
-                idx = numpy.argmin(ov)
-                act.append(idx)
-                ov[idx] = 1 * len(core)
-            act.sort()
-            return act
-
-        def get_core(active):
-            ov = numpy.zeros(shape=(len(self.integral_manager.orbitals)))
-            for i in active:
-                for j in range(len(d)):
-                    ov[j] += numpy.abs(inner(d.T[i], c.T[j], s))
-            co = []
-            for i in range(len(self.integral_manager.orbitals) - len(active)):
-                idx = numpy.argmin(ov)
-                co.append(idx)
-                ov[idx] = 1 * len(active)
-            co.sort()
-            return co
-
-        def active_to_active(active):
+            n_fr = d_orig.shape[1] - len(active_idx_c)
+            # 1. Extract the entire reference active space block from c_orig
+            c_active = c_orig[:, active_idx_c]
+
+            # 2. Compute the full overlap matrix between reference active and all d_orig orbitals
+            # Shape will be (n_active_ref, n_total_orbitals_d)
+            overlap_matrix = c_active.T @ s @ d_orig
+
+            # 3. Sum of squares along the reference axis gives the total "active character"
+            # Shape will be (n_total_orbitals_d,)
+            active_weights = numpy.sum(overlap_matrix**2, axis=0)
+
+            # 4. Sort all orbital indices of d_orig by weight in descending order
+            sorted_d_indices = numpy.argsort(active_weights)
+
+            # 5. Select the top N orbitals that match the active space best
+            chosen_active_idx = sorted_d_indices[:n_fr]
+
+            return sorted(chosen_active_idx)
+
+        def get_core(c_orig: numpy.ndarray, d_orig: numpy.ndarray, s: numpy.ndarray, active_idx_d: list[int]):
             """
-            translates active indices from canonical/the original basis to the native coeffs
+            Given the active space indices of d_orig, finds which occupied orbitals in c_orig should be frozen (core orbitals).
+
+            Parameters:
+            -----------
+            c_orig: original orbital matrix which will be frozen (typicall HF)
+            d_orig: original orbital matrix to look for active w.r.t. c_orig (i.e. native orbital matrix or CLPO matrix before active space considerations)
+            s: overlap_integrals
+            active_idx_d: The indices of the active space orbitals in d_orig.
             """
-            ov = numpy.zeros(shape=(len(self.integral_manager.orbitals)))
-            for i in active:
-                for j in range(len(d)):
-                    ov[j] += numpy.abs(inner(c.T[i], d.T[j], s))
-            act = []
-            for i in range(len(active)):
-                idx = numpy.argmax(ov)
-                act.append(idx)
-                ov[idx] = 0.0
-            act.sort()
-            return act
+            n_occ_c = d_orig.shape[1] - len(active_idx_d)
+
+            # 1. Extract the active subspace block from d_orig
+            d_active = d_orig[:, active_idx_d]
+
+            # 2. Compute the overlap between all c_orig orbitals and the d_orig active subspace
+            # Shape will be (n_total_orbitals_c, n_active_d)
+            overlap_matrix = c_orig.T @ s @ d_active
+
+            # 3. Sum of squares along the d_active axis gives the "active character" of each c_orig orbital
+            active_weights = numpy.sum(overlap_matrix**2, axis=1)
+
+            # 4. Sort the orbitals by their active weight in ASCENDING order
+            # The orbitals with the LOWEST active weight are your core (frozen) orbitals!
+            sorted_fr_indices = numpy.argsort(active_weights)
+
+            chosen_fr_idx = sorted_fr_indices[:n_occ_c]
+
+            return sorted(chosen_fr_idx)
 
         active = None
+        if not self.integral_manager.active_space_is_trivial() and core is None:
+            core = [i.idx_total for i in self.integral_manager.orbitals if i.idx is None]
         if "active" in kwargs:
             active = kwargs["active"]
             kwargs.pop("active")
-            if not len(core):
-                core = get_core(active)
+            if core is None:
+                core = get_core(c, d, s, active)
         else:
-            if not len(core):
-                if not self.integral_manager.active_space_is_trivial():
-                    active = [i.idx_total for i in self.integral_manager.orbitals if i.idx is not None]
-                    active = active_to_active(active)
-                    core = [i.idx_total for i in self.integral_manager.orbitals if i.idx is None]
-                else:
+            if active is None:
+                if core is None:
                     core = []
                     active = [i for i in range(len(self.integral_manager.orbitals))]
-            else:
-                active = get_active(core)
+                else:
+                    if isinstance(core, int):
+                        core = [core]
+                    active = get_active(c, d, s, core)
         assert len(active) + len(core) == len(self.integral_manager.orbitals)
         if "reference_orbitals" in kwargs:
             reference_orbitals = kwargs["reference_orbitals"]
@@ -831,11 +829,15 @@ def active_to_active(active):
             )
         else:
             reference_orbitals = [i.idx_total for i in self.integral_manager.reference_orbitals]
-
         to_active = [i for i in range(len(self.integral_manager.orbitals)) if i not in core]
         to_active = {active[i]: to_active[i] for i in range(len(active))}
         if len(core):
-            coeff = orthogonalize(c, d, s)
+            c_combined = numpy.zeros(shape=c.shape)
+            for i, idx in enumerate(core):
+                c_combined[:, i] = c[:, idx]
+            for act_idx in active:
+                c_combined[:, to_active[act_idx]] = d[:, act_idx]
+            coeff = orthogonalize_active_space(c_combined, s, core, [*to_active.values()])
             if inplace:
                 self.integral_manager = self.initialize_integral_manager(
                     one_body_integrals=self.integral_manager.one_body_integrals,
@@ -846,7 +848,6 @@ def active_to_active(active):
                     frozen_orbitals=core,
                     orbital_coefficients=coeff,
                     overlap_integrals=s,
-                    orbital_type="orthonormalized-{}-basis".format(self.integral_manager._basis_name),
                 )
                 return self
             else:
@@ -859,7 +860,6 @@ def active_to_active(active):
                     frozen_orbitals=core,
                     orbital_coefficients=coeff,
                     overlap_integrals=s,
-                    orbital_type="orthonormalized-{}-basis".format(self.integral_manager._basis_name),
                 )
                 parameters = copy.deepcopy(self.parameters)
                 result = QuantumChemistryBase(
diff --git a/tests/test_scipy.py b/tests/test_scipy.py
index ba64ccb4..16b08ba3 100644
--- a/tests/test_scipy.py
+++ b/tests/test_scipy.py
@@ -9,52 +9,6 @@
 samplers = select_backends.get(sampler=True)
 
 
-@pytest.mark.parametrize("simulator", simulators)
-def test_execution(simulator):
-    U = (
-        tq.gates.Rz(angle="a", target=0)
-        + tq.gates.X(target=2)
-        + tq.gates.Ry(angle="b", target=1, control=2)
-        + tq.gates.Trotterized(
-            angles=["c", "d"], generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2
-        )
-        + tq.gates.Trotterized(
-            angles=[1.0, 2.0], generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2
-        )
-        + tq.gates.ExpPauli(angle="a", paulistring="X(0)Y(1)Z(2)")
-    )
-
-    H = 1.0 * tq.paulis.X(0) + 2.0 * tq.paulis.Y(1) + 3.0 * tq.paulis.Z(2)
-    O = tq.ExpectationValue(U=U, H=H)
-
-    result = tq.optimizer_scipy.minimize(
-        objective=O, options={"maxiter": 2}, method="TNC", backend=simulator, silent=True
-    )
-
-
-@pytest.mark.parametrize("simulator", samplers)
-def test_execution_shot(simulator):
-    U = (
-        tq.gates.Rz(angle="a", target=0)
-        + tq.gates.X(target=2)
-        + tq.gates.Ry(angle="b", target=1, control=2)
-        + tq.gates.Trotterized(
-            angles=["c", "d"], generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2
-        )
-        + tq.gates.Trotterized(
-            angles=[1.0, 2.0], generators=[-0.25 * tq.paulis.Z(1), tq.paulis.X(0) + tq.paulis.Y(1)], steps=2
-        )
-        + tq.gates.ExpPauli(angle="a", paulistring="X(0)Y(1)Z(2)")
-    )
-    H = 1.0 * tq.paulis.X(0) + 2.0 * tq.paulis.Y(1) + 3.0 * tq.paulis.Z(2)
-    O = tq.ExpectationValue(U=U, H=H)
-
-    result = tq.optimizer_scipy.minimize(
-        objective=O, options={"maxiter": 2}, method="TNC", backend=simulator, samples=3, silent=True
-    )
-    assert len(result.history.energies) <= 3
-
-
 @pytest.mark.parametrize("simulator", simulators)
 def test_one_qubit_wfn(simulator):
     U = tq.gates.Trotterized(angles=["a"], steps=1, generators=[tq.paulis.Y(0)])