diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py
index 9689b042..669a3b37 100644
--- a/src/tequila/quantumchemistry/qc_base.py
+++ b/src/tequila/quantumchemistry/qc_base.py
@@ -645,6 +645,65 @@ def n_electrons(self) -> int:
         """
         return 2 * len(self.integral_manager.active_reference_orbitals)
 
+    def make_annihilation_op(self, orbital, coefficient=1.0):
+        """
+        Compute annihilation operator on orbital=orbital in qubit representation
+        """
+        assert orbital<=self.n_orbitals*2
+        aop = openfermion.ops.FermionOperator(f'{orbital}', coefficient)
+        return self.transformation(aop)
+
+    def make_creation_op(self, orbital, coefficient=1.0):
+        """
+        Compute creation operator on orbital=orbital in qubit representation
+        """
+        assert orbital<=self.n_orbitals*2
+        cop = openfermion.ops.FermionOperator(f'{orbital}^', coefficient)
+        return self.transformation(cop)
+
+    def make_number_op(self, orbital):
+        """
+        Compute number operator on orbital=orbital in qubit representation
+        """
+        num_op = self.make_creation_op(orbital) * self.make_annihilation_op(orbital)
+        return num_op
+    
+    def make_sz_op(self):
+        """
+        Compute the spin_z operator of the molecule in qubit representation
+        """
+        sz = QubitHamiltonian()
+        for i in range(0, self.n_orbitals * 2, 2):
+            one = 0.5 * self.make_creation_op(i) * self.make_annihilation_op(i)
+            two = 0.5 * self.make_creation_op(i+1) * self.make_annihilation_op(i+1)
+            sz += (one - two)
+        return sz
+
+    def make_sp_op(self):
+        """
+        Compute the spin+ operator of the molecule in qubit representation
+        """
+        sp = QubitHamiltonian()
+        for i in range(self.n_orbitals):
+            sp += self.make_creation_op(i*2) * self.make_annihilation_op(i*2 + 1)
+        return sp
+
+    def make_sm_op(self):
+        """
+        Compute the spin- operator of the molecule in qubit representation
+        """
+        sm = QubitHamiltonian()
+        for i in range(self.n_orbitals):
+            sm += self.make_creation_op(i*2 + 1) * self.make_annihilation_op(i*2)
+        return sm
+
+    def make_s2_op(self):
+        """
+        Compute the spin^2 operator of the molecule in qubit representation
+        """
+        s2_op = self.make_sm_op() * self.make_sp_op() + self.make_sz_op() * (self.make_sz_op() + 1)
+        return s2_op
+
     def make_hamiltonian(self, *args, **kwargs) -> QubitHamiltonian:
         """
         Parameters