diff --git a/.github/workflows/ci_basic.yml b/.github/workflows/ci_basic.yml
index 4b97d18..cba064d 100644
--- a/.github/workflows/ci_basic.yml
+++ b/.github/workflows/ci_basic.yml
@@ -24,7 +24,7 @@ jobs:
         run: |
           python -m pip install --upgrade pip
           pip install -r requirements.txt
-          pip install fqe
+          pip install fqe --no-deps
           pip install git+https://github.com/davibinco/TenCirChem-NG.git
       - name: Lint with flake8
         run: |
diff --git a/src/sunrise/CLPO/orbital_transformation.py b/src/sunrise/CLPO/orbital_transformation.py
index b10830c..f52a1bd 100644
--- a/src/sunrise/CLPO/orbital_transformation.py
+++ b/src/sunrise/CLPO/orbital_transformation.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 from tequila.quantumchemistry.pyscf_interface import QuantumChemistryPySCF
 import os
 import numpy
@@ -6,8 +7,6 @@
 from copy import deepcopy
 import subprocess
 from tequila.quantumchemistry.qc_base import QuantumChemistryBase
-from sunrise.molecules.hybrid_base import HybridBase
-from sunrise.molecules.fermionic_base import FermionicBase
 import numpy
 from copy import deepcopy
 from typing import Tuple
@@ -15,85 +14,9 @@
 from .binary_interface import *
 from sunrise import from_tequila
 from tequila import TequilaException
+from sunrise.molecules.utils_orbital_transformation import transform, orthogonalize
 
-def __transform(modified:QuantumChemistryBase, original:QuantumChemistryBase = None, orbital_type = 'CLPO') -> Tuple[QuantumChemistryBase, dict]:
-    '''
-    Procedure similar to what is done in use_native_orbitals but for arbitrary basis. Keeps frozen orbitals canHF
-    orthogonalized with the active modified ones
-    Returns modified molecule with the core orbitals of the original one
-    And a dictionary with the form {active_orbital_index_before:active_orbital_index_after}
-    The frozen orbitals will always be the N first on the orbital matrix
-    '''
-    def inner(a, b, s):
-        return numpy.sum(numpy.multiply(numpy.outer(a, b), s))
-    core = [i.idx_total for i in original.integral_manager.orbitals if i.idx is None]
-    assert len(original.integral_manager.orbitals) == len(modified.integral_manager.orbitals)
-    d = deepcopy(modified.integral_manager.orbital_coefficients).T
-    c = deepcopy(original.integral_manager.orbital_coefficients).T
-    s = original.integral_manager.overlap_integrals
-    n_basis = len(d)
-    ov = numpy.zeros(shape = (n_basis))
-    for i in core:
-        for j in range(n_basis):
-            ov[j] += numpy.abs(inner(c[i], d[j], s))
-    co = {}
-    for i in core:
-        idx = numpy.argmax(ov)
-        co[i] = idx
-        ov[idx] = 0
-    active = [i for i in range(n_basis) if i not in co.values()]
-    to_active =  [i for i in range(n_basis) if  i not in co.keys()]
-    to_active = {active[i] : to_active[i] for i in range(len(active))}
-    reference_orbitals = [*co.keys()]
-    i =0
-    while len(reference_orbitals) < original.parameters.total_n_electrons//2:
-        if i not in reference_orbitals:
-            reference_orbitals.append(i)
-        i += 1
-    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)
-    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
-    for i in to_active.values():
-        norm = numpy.sqrt(inner(dbar[i], dbar[i], s.T))
-        if not numpy.isclose(norm, 0):
-            dbar[i] = dbar[i] / norm
-    for j in to_active.values():
-        c[j] = dbar[j]
-    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]
-    lam_s, l_s = numpy.linalg.eigh(sprima)
-    for ei, e in enumerate(lam_s):
-        if numpy.isclose(e,0,atol=1.e-9):
-            lam_s[ei] = max(e,1.e-9) # Note: Fix to avoid inestabilities 
-    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))
-    jcoef = symm_orthog.dot(c).T
-    ref = [i.idx_total for i in original.integral_manager.reference_orbitals if i not in original.integral_manager.active_reference_orbitals]
-    ref.extend([i for i in range(n_basis) if  i not in co.keys()][:len(original.integral_manager.active_reference_orbitals)])
-    integral_manager = modified.initialize_integral_manager(one_body_integrals=original.integral_manager.one_body_integrals,
-                    two_body_integrals=original.integral_manager.two_body_integrals, constant_term=original.integral_manager.constant_term,
-                    active_orbitals= [i for i in range(n_basis) if  i not in co.keys()], frozen_orbitals=[*co.keys()], orbital_coefficients=jcoef,
-                    overlap_integrals=original.integral_manager.overlap_integrals, reference_orbitals=ref, orbital_type=orbital_type)
-    parameters = deepcopy(original.parameters)
-    if isinstance(modified, FermionicBase):
-        return FermionicBase(parameters=parameters, integral_manager=integral_manager, fermionic_backend=modified.fermionic_backend), to_active
-    elif isinstance(modified, HybridBase):
-        return HybridBase(parameters=parameters, integral_manager=integral_manager, transformation=modified.transformation, select=modified.select, two_qubit=modified.two_qubit, condense=modified.condense), to_active
-    return QuantumChemistryBase(parameters=parameters, integral_manager=integral_manager, transformation=modified.transformation), to_active
-
-def __get_MP2_occ(mol:QuantumChemistryBase):
+def __get_MP2_occ(mol:QuantumChemistryBase) -> Tuple[list[Number], list[Number]]:
     ''''
     Small helper function, given a tequila molecule returns the MP2 orbital occupation and orbital energy
     '''
@@ -103,7 +26,7 @@ def __get_MP2_occ(mol:QuantumChemistryBase):
     rdm1 = mp.MP2(hf).run().make_rdm1()
     return fr + numpy.diag(rdm1).tolist(), hf.mo_energy
 
-def generate_molden(mol:QuantumChemistryBase, filename:str = None, output_dir:str  =None, mo_occ:list = None, mo_energy:list = None, use_mp2:bool = False, option1:bool = True, use_active:bool = True):
+def generate_molden(mol:QuantumChemistryBase, filename:str = None, output_dir:str = None, mo_occ:list = None, mo_energy:list = None, use_mp2:bool = False, option1:bool = True, use_active:bool = True):
     '''
     Interface with pyscf.tools molden file generation
 
@@ -141,7 +64,7 @@ def generate_molden(mol:QuantumChemistryBase, filename:str = None, output_dir:st
     if filename is None: filename=mol.parameters.name
 
 
-    mo_coeff = mol.integral_manager.orbital_coefficients
+    mo_coeff = mol.integral_manager.orbital_coefficients.copy()
     if use_active:
         if len(mo_occ) == size_basis:
             mo_occ = [mo_occ[i] for i in active]
@@ -197,6 +120,8 @@ def generate_CLPO_molecule_edges(mol:QuantumChemistryBase, edges:list[tuple[int]
         c += f' -edges {edges}'
     call_janpa(command = c, silent = silent)
     mo_matrix = read_molden_mo_matrix(f"{output_dir}/{filename}_CLPO.molden")
+    if not use_active:
+        mo_matrix = orthogonalize(mo_matrix, mol.integral_manager.overlap_integrals)
     if rm_files:
         subprocess.call(f'rm {output_dir}/m2a.ini', shell=True)
         subprocess.call(f'rm {output_dir}/{filename}.molden', shell=True) 
@@ -205,8 +130,10 @@ def generate_CLPO_molecule_edges(mol:QuantumChemistryBase, edges:list[tuple[int]
     nmol = deepcopy(mol)
     nmol.integral_manager.orbital_coefficients = mo_matrix
     if use_active:
-        mol, to_active = __transform(original = mol, modified = nmol)
-    else: mol = nmol
+        mol, to_active = transform(original = mol, modified = nmol, orbital_type = 'CLPO')
+    else:
+        mol = nmol
+        mol.integral_manager._orbital_type = 'CLPO'
     graph = extract_clpo_graph(f"{output_dir}/graph")
     if use_active:
         ncore = len(mol.integral_manager.orbital_coefficients) - mol.n_orbitals
@@ -216,7 +143,7 @@ def generate_CLPO_molecule_edges(mol:QuantumChemistryBase, edges:list[tuple[int]
         subprocess.call(f'rm {output_dir}/graph', shell=True)
     return mol,graph
 
-def generate_HAO_molecule(mol:QuantumChemistryBase, output_dir:str = None, thres:Number = 1.e-9, silent:bool = True, use_active:bool = True, rm_files:bool = True,**kwargs)->QuantumChemistryBase:
+def generate_HAO_molecule(mol:QuantumChemistryBase, output_dir:str = None, thres:Number = 1.e-9, silent:bool = True, use_active:bool = True, rm_files:bool = True,**kwargs) -> QuantumChemistryBase:
     '''
     Temporal function for generating a molecule with Hybrid Atomic Orbitals via janpa (10.1002/qua.25798) until integrated in Sunrise molecules
     
@@ -245,6 +172,8 @@ def generate_HAO_molecule(mol:QuantumChemistryBase, output_dir:str = None, thres
     call_molden2molden(command = f'-NormalizeBF -cart2pure  -i {filename}.molden -o {filename}.molden', silent = silent, output_dir = output_dir)
     call_janpa(command=f'-i {filename}.molden -AHO_Molden_File {filename}_HAO.molden -HybrOptOccConvThresh {thres}', silent = silent, output_dir = output_dir)
     mo_matrix = read_molden_mo_matrix(f"{output_dir}/{filename}_HAO.molden")
+    if not use_active:
+        mo_matrix = orthogonalize(mo_matrix, mol.integral_manager.overlap_integrals)
     if rm_files:
         subprocess.call(f'rm {output_dir}/m2a.ini', shell=True)
         subprocess.call(f'rm {output_dir}/{filename}.molden', shell=True)
@@ -254,8 +183,10 @@ def generate_HAO_molecule(mol:QuantumChemistryBase, output_dir:str = None, thres
     nmol = deepcopy(mol)
     nmol.integral_manager.orbital_coefficients = mo_matrix
     if use_active:
-        mol, to_active = __transform(original = mol, modified = nmol, orbital_type = 'HAO')
-    else: mol = nmol
+        mol, to_active = transform(original = mol, modified = nmol, orbital_type = 'HAO')
+    else:
+        mol = nmol
+        mol.integral_manager._orbital_type = "HAO"
     if ret2act:
         return mol, to_active
     return mol
diff --git a/src/sunrise/MCVBT/GNM.py b/src/sunrise/MCVBT/GNM.py
index 1ee99d8..113bdb4 100644
--- a/src/sunrise/MCVBT/GNM.py
+++ b/src/sunrise/MCVBT/GNM.py
@@ -8,7 +8,10 @@
 
 import sunrise as sn
 from sunrise.MCVBT.Big_Exp import BigExpVal
-from sunrise.expval.fqe_expval import FQEBraKet
+try:
+    from sunrise.expval.fqe_expval import FQEBraKet
+except ImportError:
+    pass
 
 from typing import Dict
 import csv
diff --git a/src/sunrise/expval/fqe_circuit_sim.py b/src/sunrise/expval/fqe_circuit_sim.py
index dcac528..bfa6184 100644
--- a/src/sunrise/expval/fqe_circuit_sim.py
+++ b/src/sunrise/expval/fqe_circuit_sim.py
@@ -2,7 +2,10 @@
 from tequila import Variable,Objective,simulate,QubitWaveFunction,TequilaWarning,BitNumbering
 from tequila.objective.objective import Variables
 from tequila import  BitString, BitStringLSB
-from sunrise.expval.fqe_expval import FQEBraKet
+try:
+    from sunrise.expval.fqe_expval import FQEBraKet
+except ImportError:
+    pass
 from sunrise.fermionic_operations import FCircuit
 from tequila import Molecule
 import numpy as np
diff --git a/src/sunrise/molecules/fermionic_base/fer_base.py b/src/sunrise/molecules/fermionic_base/fer_base.py
index 6c3d85f..e8ecf53 100644
--- a/src/sunrise/molecules/fermionic_base/fer_base.py
+++ b/src/sunrise/molecules/fermionic_base/fer_base.py
@@ -438,84 +438,10 @@ def use_native_orbitals(self, inplace=False, core: list = None, *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()
-
-        def inner(a, b, s):
-            return numpy.sum(numpy.multiply(numpy.outer(a, b), s))
-
-        def orthogonalize(c, d, s):
-            """
-            :return: orthogonalized orbitals with core HF orbitals and active Orthongonalized Native orbitals.
-            """
-            ### 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
+        from sunrise.molecules.utils_orbital_transformation import get_active, get_core, orthogonalize_active_space
+        c = self.integral_manager.orbital_coefficients.copy()
+        s = self.integral_manager.overlap_integrals.copy()
+        d = self.integral_manager.get_orthonormalized_orbital_coefficients().copy()
 
         active = None
         if not self.integral_manager.active_space_is_trivial() and core is None:
@@ -524,7 +450,7 @@ def get_core(active):
             active = kwargs["active"]
             kwargs.pop("active")
             if core is None:
-                core = get_core(active)
+                core = get_core(c, d, s, active)
         else:
             if active is None:
                 if core is None:
@@ -533,19 +459,30 @@ def get_core(active):
                 else:
                     if isinstance(core, int):
                         core = [core]
-                    active = get_active(core)
+                    active = get_active(c, d, s, [i.idx_total for i in self.integral_manager.active_orbitals])
         assert len(active) + len(core) == len(self.integral_manager.orbitals)
+        if "reference_orbitals" in kwargs:
+            reference_orbitals = kwargs["reference_orbitals"]
+            kwargs.pop()
+            assert len(reference_orbitals) == len(self.parameters.total_n_electrons)//2,f'Number of  provided reference_orbitals incorrect. Expected {self.parameters.total_n_electrons//2}, received {len(reference_orbitals)}'
+        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,
                     two_body_integrals=self.integral_manager.two_body_integrals,
                     constant_term=self.integral_manager.constant_term,
                     active_orbitals=[*to_active.values()],
-                    reference_orbitals=[i.idx_total for i in self.integral_manager.reference_orbitals],
+                    reference_orbitals=reference_orbitals,
                     frozen_orbitals=core,
                     orbital_coefficients=coeff,
                     overlap_integrals=s,
diff --git a/src/sunrise/molecules/hybrid_base/HybridBase.py b/src/sunrise/molecules/hybrid_base/HybridBase.py
index 5e2acdd..cd3d731 100644
--- a/src/sunrise/molecules/hybrid_base/HybridBase.py
+++ b/src/sunrise/molecules/hybrid_base/HybridBase.py
@@ -20,7 +20,7 @@
 from openfermion import FermionOperator
 import copy
 from sunrise.hybridization.hybridization import Graph
-from typing import Union,Optional,List
+from typing import Union, Optional, List
 class HybridBase(qc_base):
     def __init__(self, parameters: ParametersQC,select: typing.Union[str,dict]={},transformation: typing.Union[str, typing.Callable] = None, active_orbitals: list = None,
                  frozen_orbitals: list = None, orbital_type: str = None,reference_orbitals: list = None, orbitals: list = None, *args, **kwargs):
@@ -257,125 +257,44 @@ def use_native_orbitals(self, inplace=False, core: list = None, *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()
-
-        def inner(a, b, s):
-            return numpy.sum(numpy.multiply(numpy.outer(a, b), s))
-
-        def orthogonalize(c, d, s):
-            '''
-            :return: orthogonalized orbitals with core HF orbitals and active Orthongonalized Native orbitals.
-            '''
-            ### 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):
-            '''
-            translates active indices from canonical/the original basis to the native coeffs
-            '''
-            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.
-            act.sort()
-            return act
+        from sunrise.molecules.utils_orbital_transformation import get_active, get_core, orthogonalize_active_space
+        c = self.integral_manager.orbital_coefficients.copy()
+        s = self.integral_manager.overlap_integrals.copy()
+        d = self.integral_manager.get_orthonormalized_orbital_coefficients().copy()
 
         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 core is None:
-                core = get_core(active)
+                core = get_core(c, d, s, active)
         else:
-            if core is None:
-                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:
-                if isinstance(core, int):
-                    core = [core]
-                active = get_active(core)
+                else:
+                    if isinstance(core, int):
+                        core = [core]
+                    active = get_active(c, d, s, [i.idx_total for i in self.integral_manager.active_orbitals])
         assert len(active) + len(core) == len(self.integral_manager.orbitals)
+        if "reference_orbitals" in kwargs:
+            reference_orbitals = kwargs["reference_orbitals"]
+            kwargs.pop()
+            assert len(reference_orbitals) == len(self.parameters.total_n_electrons)//2,f'Number of  provided reference_orbitals incorrect. Expected {self.parameters.total_n_electrons//2}, received {len(reference_orbitals)}'
+        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 not all([i == to_active[i] for i in to_active]) and len(self.BOS_MO) and len(self.FER_MO):
                 print("Orbital may be reordered, please double check F/B selection")
             if len(active) == len(self.select):
@@ -389,8 +308,8 @@ def active_to_active(active):
                     two_body_integrals=self.integral_manager.two_body_integrals,
                     constant_term=self.integral_manager.constant_term,
                     active_orbitals=[*to_active.values()],
-                    reference_orbitals=[i.idx_total for i in self.integral_manager.reference_orbitals]
-                    , frozen_orbitals=core, orbital_coefficients=coeff, overlap_integrals=s,
+                    reference_orbitals=reference_orbitals,
+                    frozen_orbitals=core, orbital_coefficients=coeff, overlap_integrals=s,
                     orbital_type="orthonormalized-{}-basis".format(self.integral_manager._basis_name),
                     )
                 self.update_select(new_select)
@@ -399,10 +318,10 @@ def active_to_active(active):
                 integral_manager = self.initialize_integral_manager(
                     one_body_integrals=self.integral_manager.one_body_integrals,
                     two_body_integrals=self.integral_manager.two_body_integrals,
-                    constant_term=self.integral_manager.constant_term
-                    , active_orbitals=[*to_active.values()],
-                    reference_orbitals=[i.idx_total for i in self.integral_manager.reference_orbitals]
-                    , frozen_orbitals=core, orbital_coefficients=coeff, overlap_integrals=s,
+                    constant_term=self.integral_manager.constant_term,
+                    active_orbitals=[*to_active.values()],
+                    reference_orbitals=reference_orbitals,
+                    frozen_orbitals=core, orbital_coefficients=coeff, overlap_integrals=s,
                     orbital_type="orthonormalized-{}-basis".format(self.integral_manager._basis_name),
                     )
                 parameters = copy.deepcopy(self.parameters)
diff --git a/src/sunrise/molecules/utils_orbital_transformation.py b/src/sunrise/molecules/utils_orbital_transformation.py
new file mode 100644
index 0000000..6eacf2e
--- /dev/null
+++ b/src/sunrise/molecules/utils_orbital_transformation.py
@@ -0,0 +1,159 @@
+import numpy
+from typing import Tuple
+from copy import deepcopy
+from tequila.quantumchemistry.qc_base import QuantumChemistryBase
+from sunrise.molecules.hybrid_base import HybridBase
+from sunrise.molecules.fermionic_base import FermionicBase
+
+def orthogonalize(c:numpy.ndarray, s:numpy.ndarray) -> numpy.ndarray:
+    """
+    Symmetrically orthogonalize orbital coefficients.
+    c: (basis_functions, orbitals)
+    s: (basis_functions, basis_functions)
+    """
+    # 1. Compute the overlap of the current MOs: S' = C^T * S * C
+    # This replaces your entire nested loop and inner() function.
+    sprima = c.T @ s @ c
+    # 2. Diagonalize S'
+    lam_s, l_s = numpy.linalg.eigh(sprima)
+
+    # Optional but recommended: Clip tiny negative eigenvalues due to numerical noise
+    lam_s = numpy.maximum(lam_s, 1e-14)
+
+    # 3. Construct (S')^{-1/2}
+    # This is much faster/cleaner than inverting a full matrix
+    lam_sqrt_inv = numpy.diag(1.0 / numpy.sqrt(lam_s))
+    symm_orthog = l_s @ lam_sqrt_inv @ l_s.T
+
+    # 4. Transform coefficients: C_new = C * (S')^{-1/2}
+    jcoef = c @ symm_orthog
+
+    return jcoef
+
+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 get_active(c_orig:numpy.ndarray, d_orig:numpy.ndarray, s:numpy.ndarray, active_idx_c:list[int]) -> list[int]:
+    """
+    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
+    """
+    # 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)[::-1]
+
+    # 5. Select the top N orbitals that match the active space best
+    chosen_active_idx = sorted_d_indices[:len(active_idx_c)]
+
+    return sorted(chosen_active_idx)
+
+def get_core(c_orig:numpy.ndarray, d_orig:numpy.ndarray, s:numpy.ndarray, active_idx_d:list[int]):
+    """
+    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.
+    """
+    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)
+
+def transform(modified:QuantumChemistryBase, original:QuantumChemistryBase, orbital_type:str = None) -> Tuple[QuantumChemistryBase, dict]:
+    '''
+    Procedure similar to what is done in use_native_orbitals but for arbitrary basis keeped insied modified. Keeps frozen orbitals canHF
+    orthogonalized with the active modified ones
+    Returns modified molecule with the core orbitals of the original one
+    And a dictionary with the form {active_orbital_index_before:active_orbital_index_after}
+    The frozen orbitals will always be the N first on the orbital matrix
+    '''
+    core = [i.idx_total for i in original.integral_manager.orbitals if i.idx is None]
+    assert len(original.integral_manager.orbitals) == len(modified.integral_manager.orbitals)
+    c_orig = original.integral_manager.orbital_coefficients.copy()
+    d_orig = modified.integral_manager.orbital_coefficients.copy()
+    s = original.integral_manager.overlap_integrals.copy()
+    n_basis = c_orig.shape[0]
+    active = get_active(c_orig, d_orig, s, [i.idx_total for i in original.integral_manager.orbitals if i.idx is not None])
+    to_active =  [i for i in range(n_basis) if  i not in core]
+    to_active = {active[i] : to_active[i] for i in range(len(active))}
+    reference_orbitals = core.copy()
+    i =0
+    while len(reference_orbitals) < original.parameters.total_n_electrons//2:
+        if i not in reference_orbitals:
+            reference_orbitals.append(i)
+        i += 1
+    n_core = len(core)
+    n_active = len(active)
+    c_combined = numpy.zeros((n_basis, n_core + n_active))
+    for i,idx in enumerate(core):
+        c_combined[:, i] = c_orig[:, idx]
+    for act_idx in active:
+        c_combined[:, to_active[act_idx]] = d_orig[:, act_idx]
+    jcoef = orthogonalize_active_space(c_combined, s,core, [*to_active.values()])
+    ref = [i.idx_total for i in original.integral_manager.reference_orbitals if i not in original.integral_manager.active_reference_orbitals]
+    ref.extend([i for i in range(n_basis) if  i not in core][:len(original.integral_manager.active_reference_orbitals)])
+    integral_manager = modified.initialize_integral_manager(one_body_integrals=original.integral_manager.one_body_integrals,
+                    two_body_integrals=original.integral_manager.two_body_integrals, constant_term=original.integral_manager.constant_term,
+                    active_orbitals= [i for i in range(n_basis) if  i not in core], frozen_orbitals=core, orbital_coefficients=jcoef,
+                    overlap_integrals=original.integral_manager.overlap_integrals, reference_orbitals=ref, orbital_type=orbital_type)
+    parameters = deepcopy(original.parameters)
+    if isinstance(modified, FermionicBase):
+        return FermionicBase(parameters=parameters, integral_manager=integral_manager, fermionic_backend=modified.fermionic_backend), to_active
+    elif isinstance(modified, HybridBase):
+        return HybridBase(parameters=parameters, integral_manager=integral_manager, transformation=modified.transformation, select=modified.select, two_qubit=modified.two_qubit, condense=modified.condense), to_active
+    return QuantumChemistryBase(parameters=parameters, integral_manager=integral_manager, transformation=modified.transformation), to_active