diff --git a/quantum_systems/custom_system.py b/quantum_systems/custom_system.py
index 9efb441..dc31c80 100644
--- a/quantum_systems/custom_system.py
+++ b/quantum_systems/custom_system.py
@@ -111,6 +111,8 @@ def construct_pyscf_system_rhf(
     verbose=False,
     charge=0,
     cart=False,
+    nfrozen=0,
+    conv_tol=1e-10,
     **kwargs,
 ):
     """Convenience function setting up a closed-shell atom or a molecule from
@@ -187,6 +189,7 @@ def construct_pyscf_system_rhf(
     l = mol.nao
 
     hf = pyscf.scf.RHF(mol)
+    hf.conv_tol=conv_tol
     hf_energy = hf.kernel()
 
     if not hf.converged:
@@ -215,7 +218,7 @@ def construct_pyscf_system_rhf(
     bs.dipole_moment = dipole_integrals
     bs.change_module(np=np)
 
-    system = SpatialOrbitalSystem(n, bs)
+    system = SpatialOrbitalSystem(n, bs, nfrozen=nfrozen)
     system.change_basis(C)
 
     return (
diff --git a/quantum_systems/general_orbital_system.py b/quantum_systems/general_orbital_system.py
index 4f2f25e..afddbcd 100644
--- a/quantum_systems/general_orbital_system.py
+++ b/quantum_systems/general_orbital_system.py
@@ -55,7 +55,7 @@ def compute_reference_energy(self):
             The reference energy.
         """
 
-        o, v = self.o, self.v
+        o, v = slice(0, self.n), self.v
 
         return self.np.trace(self.h[o, o]) + 0.5 * self.np.trace(
             self.np.trace(self.u[o, o, o, o], axis1=1, axis2=3)
@@ -91,7 +91,7 @@ def construct_fock_matrix(self, h, u, f=None):
         """
 
         np = self.np
-        o, v = (self.o, self.v)
+        o, v = (slice(0, self.n), self.v)
 
         if f is None:
             f = np.zeros_like(h)
diff --git a/quantum_systems/spatial_orbital_system.py b/quantum_systems/spatial_orbital_system.py
index 0ab0da1..1d90048 100644
--- a/quantum_systems/spatial_orbital_system.py
+++ b/quantum_systems/spatial_orbital_system.py
@@ -109,7 +109,7 @@ def compute_reference_energy(self):
             The reference energy.
         """
 
-        o, v = self.o, self.v
+        o, v = slice(0, self.n), self.v
 
         return (
             2 * self.np.trace(self.h[o, o])
@@ -147,7 +147,7 @@ def construct_fock_matrix(self, h, u, f=None):
             The filled Fock matrix.
         """
         np = self.np
-        o, v = (self.o, self.v)
+        o, v = (slice(0, self.n), self.v)
 
         if f is None:
             f = np.zeros_like(h)
diff --git a/quantum_systems/system.py b/quantum_systems/system.py
index 0856ad8..d5df40e 100644
--- a/quantum_systems/system.py
+++ b/quantum_systems/system.py
@@ -16,17 +16,17 @@ class QuantumSystem(metaclass=abc.ABCMeta):
         states.
     """
 
-    def __init__(self, n, basis_set):
+    def __init__(self, n, basis_set, nfrozen=0):
         self._basis_set = basis_set
 
         assert n <= self._basis_set.l
 
         self.np = self._basis_set.np
-        self.set_system_size(n, self._basis_set.l)
+        self.set_system_size(n, self._basis_set.l, nfrozen)
 
         self._time_evolution_operator = None
 
-    def set_system_size(self, n, l):
+    def set_system_size(self, n, l, nfrozen):
         """Function setting the system size. Note that ``l`` should
         correspond to the length of each axis of the matrix elements.
 
@@ -43,8 +43,9 @@ def set_system_size(self, n, l):
         self.n = n
         self.l = l
         self.m = self.l - self.n
+        self.nfrozen = nfrozen
 
-        self.o = slice(0, self.n)
+        self.o = slice(nfrozen, self.n)
         self.v = slice(self.n, self.l)
 
     @abc.abstractmethod
@@ -65,7 +66,7 @@ def change_module(self, np):
 
     def change_basis(self, C, C_tilde=None):
         self._basis_set.change_basis(C, C_tilde)
-        self.set_system_size(self.n, self._basis_set.l)
+        self.set_system_size(self.n, self._basis_set.l, self.nfrozen)
 
     @abc.abstractmethod
     def change_to_hf_basis(self, *args, **kwargs):