Repository: PennyLaneAI/pennylane
Pull request #5638 (2024-05-03 13:35:28)

New branch for qchem features

**Context:**
A new branch for quantum chemistry features

**Description of the Change:**

**Benefits:**
Helps add small changes to this branch before merging to master

**Possible Drawbacks:**

**Related GitHub Issues:**


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5638?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.69%. Comparing base [(`b556393`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/b556393f129deb5f8c449c178428c0aadd51ff2c?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`283c7af`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5638?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff            @@
##           master    #5638    +/-   ##
========================================
  Coverage   99.69%   99.69%            
========================================
  Files         412      412            
  Lines       38599    38308   -291     
========================================
- Hits        38480    38190   -290     
+ Misses        119      118     -1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5638?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/code/qml_qchem.rst b/doc/code/qml_qchem.rst
index 68103da28..3e93340d6 100644
--- a/doc/code/qml_qchem.rst
+++ b/doc/code/qml_qchem.rst
@@ -59,7 +59,8 @@ We then construct the Hamiltonian.
 
     args = [geometry, alpha, coeff] # initial values of the differentiable parameters
 
-    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, geometry, alpha=alpha, coeff=coeff, args=args)
+    molecule = qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
+    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(molecule, args=args)
 
 >>> print(hamiltonian)
   (-0.35968235922631075) [I0]
@@ -95,7 +96,7 @@ molecular geometry optimization with PennyLane is provided in this
         def circuit(*args):
             qml.BasisState(hf_state, wires=[0, 1, 2, 3])
             qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])
-            return qml.expval(qml.qchem.molecular_hamiltonian(mol.symbols, mol.coordinates, alpha=mol.alpha, coeff=mol.coeff, args=args[1:])[0])
+            return qml.expval(qml.qchem.molecular_hamiltonian(mol, args=args[1:])[0])
         return circuit
 
 Now that the circuit is defined, we can create a geometry and parameter optimization loop. For
@@ -170,11 +171,12 @@ integral can be differentiated with respect to the basis set parameters as follo
 OpenFermion-PySCF backend
 -------------------------
 
-The :func:`~.molecular_hamiltonian` function can be also used to construct the molecular Hamiltonian
-with a non-differentiable backend that uses the
-`OpenFermion-PySCF <https://github.com/quantumlib/OpenFermion-PySCF>`_ plugin interfaced with the
+The :func:`~.molecular_hamiltonian` function can also be used to construct the molecular Hamiltonian
+with non-differentiable backends that use the
+`OpenFermion-PySCF <https://github.com/quantumlib/OpenFermion-PySCF>`_ plugin or the
 electronic structure package `PySCF <https://github.com/sunqm/pyscf>`_. The non-differentiable
-backend can be selected by setting ``method='pyscf'`` in :func:`~.molecular_hamiltonian`:
+backends can be selected by setting ``method='openfermion'`` or ``method='pyscf'``
+in ``molecular_hamiltonian``:
 
 .. code-block:: python
 
@@ -183,18 +185,13 @@ backend can be selected by setting ``method='pyscf'`` in :func:`~.molecular_hami
 
     symbols = ["H", "H"]
     geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
-    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(
-        symbols,
-        geometry,
-        charge=0,
-        mult=1,
-        basis='sto-3g',
-        method='pyscf'
-    )
-
-The non-differentiable backend requires the ``OpenFermion-PySCF`` plugin to be installed by the user
-with
+    molecule = qml.qchem.Molecule(symbols, geometry, charge=0, mult=1, basis_name='sto-3g')
+    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(molecule, method='pyscf')
+
+The non-differentiable backends require either ``OpenFermion-PySCF`` or ``PySCF`` to be installed by
+the user with
 
 .. code-block:: bash
 
     pip install openfermionpyscf
+    pip install pyscf
\ No newline at end of file
diff --git a/doc/introduction/chemistry.rst b/doc/introduction/chemistry.rst
index c37772a81..43ae07fd4 100644
--- a/doc/introduction/chemistry.rst
+++ b/doc/introduction/chemistry.rst
@@ -27,7 +27,8 @@ to generate the electronic Hamiltonian in a single call. For example,
 
     symbols = ["H", "H"]
     geometry = np.array([[0., 0., -0.66140414], [0., 0., 0.66140414]])
-    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, geometry)
+    molecule = qml.qchem.Molecule(symbols, geometry)
+    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(molecule)
 
 where:
 
@@ -36,26 +37,28 @@ where:
 * ``qubits`` is the number of qubits needed to perform the quantum simulation.
 
 The :func:`~.molecular_hamiltonian` function can also be used to construct the molecular Hamiltonian
-with an external backend that uses the
-`OpenFermion-PySCF <https://github.com/quantumlib/OpenFermion-PySCF>`_ plugin interfaced with the
+with external backends that use the
+`OpenFermion-PySCF <https://github.com/quantumlib/OpenFermion-PySCF>`_ plugin or the
 electronic structure package `PySCF <https://github.com/pyscf/pyscf>`_, which requires separate
-installation. This backend is non-differentiable and can be selected by setting
-``method='pyscf'`` in :func:`~.molecular_hamiltonian`. 
+installation. These backends are non-differentiable and can be selected by setting
+``method='openfermion'`` or ``method='pyscf'`` in ``molecular_hamiltonian``.
 
 Furthermore, the net charge,
 the `spin multiplicity <https://en.wikipedia.org/wiki/Multiplicity_(chemistry)>`_, the
-`atomic basis functions <https://www.basissetexchange.org/>`_ and the active space can also be
-specified for each backend.
+`atomic basis functions <https://www.basissetexchange.org/>`_, the mapping method and the active
+space can also be specified for each backend.
 
 .. code-block:: python
 
-    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(
+    molecule = qml.qchem.Molecule(
         symbols,
         geometry,
         charge=0,
         mult=1,
-        basis='sto-3g',
-        method='pyscf',
+        basis_name='sto-3g')
+    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(
+        molecule,
+        mapping='jordan_wigner',
         active_electrons=2,
         active_orbitals=2
     )
@@ -67,7 +70,7 @@ If the electronic Hamiltonian is built independently using
 `OpenFermion <https://github.com/quantumlib/OpenFermion>`_ tools, it can be readily converted 
 to a PennyLane observable using the :func:`~.pennylane.import_operator` function. There is also 
 capability to import wavefunctions (states) that have been pre-computed by traditional quantum chemistry methods
-from `PySCF <https://github.com/pyscf/pyscf>`_, which could be used to for example to provide a better
+from `PySCF <https://github.com/pyscf/pyscf>`_, which could be used to for example provide a better
 starting point to a quantum algorithm. State import can be accomplished using the :func:`~pennylane.qchem.import_state` 
 utility function.
 
@@ -100,7 +103,8 @@ expectation value of a Hamiltonian can be calculated using ``qml.expval``:
 
     symbols = ["H", "H"]
     geometry = np.array([[0., 0., -0.66140414], [0., 0., 0.66140414]])
-    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, geometry)
+    molecule = qml.qchem.Molecule(symbols, geometry)
+    hamiltonian, qubits = qml.qchem.molecular_hamiltonian(molecule)
 
     @qml.qnode(dev)
     def circuit(params):
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 0e4fc1531..6c9c948f4 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -111,6 +111,22 @@
   `par_info`, `obs_sharing_wires`, and `obs_sharing_wires_id` are now public attributes.
   [(#5696)](https://github.com/PennyLaneAI/pennylane/pull/5696)
 
+* The `qml.qchem.Molecule` object is now the central object used by all qchem functions.
+  [(#5571)](https://github.com/PennyLaneAI/pennylane/pull/5571)
+
+* The `qml.qchem.Molecule` class now supports Angstrom as a unit.
+  [(#5694)](https://github.com/PennyLaneAI/pennylane/pull/5694)
+
+* The `qml.qchem.Molecule` class now supports open-shell systems.
+  [(#5655)](https://github.com/PennyLaneAI/pennylane/pull/5655)
+
+* The `qml.qchem.molecular_hamiltonian` function now supports parity and Bravyi-Kitaev mappings.
+  [(#5657)](https://github.com/PennyLaneAI/pennylane/pull/5657/)
+
+* The qchem docs are updated with the new qchem improvements.
+  [(#5758)](https://github.com/PennyLaneAI/pennylane/pull/5758/)
+  [(#5638)](https://github.com/PennyLaneAI/pennylane/pull/5638/)
+
 <h4>Community contributions 🥳</h4>
 
 * Implemented kwargs (`check_interface`, `check_trainability`, `rtol` and `atol`) support in `qml.equal` for the operators `Pow`, `Adjoint`, `Exp`, and `SProd`.
@@ -219,6 +235,7 @@ Gabriel Bottrill,
 Astral Cai,
 Ahmed Darwish,
 Isaac De Vlugt,
+Diksha Dhawan,
 Pietropaolo Frisoni,
 Emiliano Godinez,
 David Ittah,
diff --git a/pennylane/qchem/hamiltonian.py b/pennylane/qchem/hamiltonian.py
index 629bc229a..f0d6696a4 100644
--- a/pennylane/qchem/hamiltonian.py
+++ b/pennylane/qchem/hamiltonian.py
@@ -181,7 +181,7 @@ def fermionic_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None):
     return _fermionic_hamiltonian
 
 
-def diff_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None):
+def diff_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None, mapping="jordan_wigner"):
     r"""Return a function that computes the qubit Hamiltonian.
 
     Args:
@@ -189,6 +189,8 @@ def diff_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None):
         cutoff (float): cutoff value for discarding the negligible electronic integrals
         core (list[int]): indices of the core orbitals
         active (list[int]): indices of the active orbitals
+        mapping (str): Specifies the fermion-to-qubit mapping. Input values can
+            be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``.
 
     Returns:
         function: function that computes the qubit hamiltonian
@@ -222,6 +224,6 @@ def diff_hamiltonian(mol, cutoff=1.0e-12, core=None, active=None):
 
         h_ferm = fermionic_hamiltonian(mol, cutoff, core, active)(*args)
 
-        return qubit_observable(h_ferm)
+        return qubit_observable(h_ferm, mapping=mapping)
 
     return _molecular_hamiltonian
diff --git a/pennylane/qchem/hartree_fock.py b/pennylane/qchem/hartree_fock.py
index 381183372..69750cc8d 100644
--- a/pennylane/qchem/hartree_fock.py
+++ b/pennylane/qchem/hartree_fock.py
@@ -116,6 +116,11 @@ def scf(mol, n_steps=50, tol=1e-8):
             tuple(array[float]): eigenvalues of the Fock matrix, molecular orbital coefficients,
             Fock matrix, core matrix
         """
+        if mol.n_electrons % 2 == 1 or mol.mult != 1:
+            raise ValueError(
+                "Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule."
+            )
+
         basis_functions = mol.basis_set
         charges = mol.nuclear_charges
         r = mol.coordinates
diff --git a/pennylane/qchem/molecule.py b/pennylane/qchem/molecule.py
index 94a2fb101..c006e4bd0 100644
--- a/pennylane/qchem/molecule.py
+++ b/pennylane/qchem/molecule.py
@@ -27,6 +27,9 @@ from .basis_data import atomic_numbers
 from .basis_set import BasisFunction, mol_basis_data
 from .integrals import contracted_norm, primitive_norm
 
+# Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0)
+bohr_angs = 0.529177210903
+
 
 class Molecule:
     r"""Create a molecule object that stores molecular information and default basis set parameters.
@@ -41,16 +44,17 @@ class Molecule:
             where ``N`` is the number of atoms.
         charge (int): net charge of the molecule
         mult (int): Spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1` for
-            :math:`N_\mathrm{unpaired}` unpaired electrons occupying the HF orbitals. Currently,
-            openshell systems are not supported; ``mult`` must be equal to :math:`1`.
+            :math:`N_\mathrm{unpaired}` unpaired electrons occupying the HF orbitals.
         basis_name (str): Atomic basis set used to represent the molecular orbitals. Currently, the
-            only supported basis sets are 'STO-3G', '6-31G', '6-311G' and 'CC-PVDZ'.
+            only supported basis sets are ``STO-3G``, ``6-31G``, ``6-311G`` and ``CC-PVDZ``. Other
+            basis sets can be loaded from the basis-set-exchange library using ``load_data``.
         load_data (bool): flag to load data from the basis-set-exchange library
         l (tuple[int]): angular momentum quantum numbers of the basis function
         alpha (array[float]): exponents of the primitive Gaussian functions
         coeff (array[float]): coefficients of the contracted Gaussian functions
         r (array[float]): positions of the Gaussian functions
         normalize (bool): if True, the basis functions get normalized
+        unit (str): unit of atomic coordinates. Available options are ``unit="bohr"`` and ``unit="angstrom"``.
 
     **Example**
 
@@ -69,11 +73,13 @@ class Molecule:
         charge=0,
         mult=1,
         basis_name="sto-3g",
+        name="molecule",
         load_data=False,
         l=None,
         alpha=None,
         coeff=None,
         normalize=True,
+        unit="bohr",
     ):
         if (
             basis_name.lower()
@@ -96,22 +102,27 @@ class Molecule:
 
         self.symbols = symbols
         self.coordinates = coordinates
+        self.unit = unit.strip().lower()
         self.charge = charge
         self.mult = mult
         self.basis_name = basis_name.lower()
-
+        self.name = name
+        self.load_data = load_data
         self.n_basis, self.basis_data = mol_basis_data(self.basis_name, self.symbols, load_data)
 
+        if self.unit not in ("angstrom", "bohr"):
+            raise ValueError(
+                f"The provided unit '{unit}' is not supported. "
+                f"Please set 'unit' to 'bohr' or 'angstrom'."
+            )
+
+        if self.unit == "angstrom":
+            self.coordinates = self.coordinates / bohr_angs
+
         self.nuclear_charges = [atomic_numbers[s] for s in self.symbols]
 
         self.n_electrons = sum(self.nuclear_charges) - self.charge
 
-        if self.n_electrons % 2 == 1 or self.mult != 1:
-            raise ValueError(
-                "Openshell systems are not supported. Change the charge or spin "
-                "multiplicity of the molecule."
-            )
-
         if l is None:
             l = [i[0] for i in self.basis_data]
 
diff --git a/pennylane/qchem/observable_hf.py b/pennylane/qchem/observable_hf.py
index ffc5326d2..bc18dc182 100644
--- a/pennylane/qchem/observable_hf.py
+++ b/pennylane/qchem/observable_hf.py
@@ -98,13 +98,14 @@ def fermionic_observable(constant, one=None, two=None, cutoff=1.0e-12):
     return sentence
 
 
-def qubit_observable(o_ferm, cutoff=1.0e-12):
+def qubit_observable(o_ferm, cutoff=1.0e-12, mapping="jordan_wigner"):
     r"""Convert a fermionic observable to a PennyLane qubit observable.
 
     Args:
         o_ferm (Union[FermiWord, FermiSentence]): fermionic operator
         cutoff (float): cutoff value for discarding the negligible terms
-
+        mapping (str): Specifies the fermion-to-qubit mapping. Input values can
+            be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``.
     Returns:
         Operator: Simplified PennyLane Hamiltonian
 
@@ -132,7 +133,15 @@ def qubit_observable(o_ferm, cutoff=1.0e-12):
     + ((0.775+0j)) [Y1 Y2]
     + ((0.775+0j)) [X1 X2]
     """
-    h = qml.jordan_wigner(o_ferm, ps=True, tol=cutoff)
+    if mapping == "jordan_wigner":
+        h = qml.jordan_wigner(o_ferm, ps=True, tol=cutoff)
+    elif mapping == "parity":
+        qubits = len(o_ferm.wires)
+        h = qml.parity_transform(o_ferm, qubits, ps=True, tol=cutoff)
+    elif mapping == "bravyi_kitaev":
+        qubits = len(o_ferm.wires)
+        h = qml.bravyi_kitaev(o_ferm, qubits, ps=True, tol=cutoff)
+
     h.simplify(tol=cutoff)
 
     if active_new_opmath():
diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py
index 1f049ef76..6373668ee 100644
--- a/pennylane/qchem/openfermion_obs.py
+++ b/pennylane/qchem/openfermion_obs.py
@@ -16,12 +16,16 @@
 # pylint: disable=too-many-arguments, too-few-public-methods, too-many-branches, unused-variable
 # pylint: disable=consider-using-generator, protected-access
 import os
+from functools import singledispatch
 
 import numpy as np
 
 import pennylane as qml
 from pennylane.operation import active_new_opmath
 from pennylane.pauli.utils import simplify
+from pennylane.qchem.molecule import Molecule
+
+from .basis_data import atomic_numbers
 
 # Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0)
 bohr_angs = 0.529177210903
@@ -592,19 +596,6 @@ def dipole_of(
     """
     openfermion, _ = _import_of()
 
-    atomic_numbers = {
-        "H": 1,
-        "He": 2,
-        "Li": 3,
-        "Be": 4,
-        "B": 5,
-        "C": 6,
-        "N": 7,
-        "O": 8,
-        "F": 9,
-        "Ne": 10,
-    }
-
     if mult != 1:
         raise ValueError(
             f"Currently, this functionality is constrained to Hartree-Fock states with spin multiplicity = 1;"
@@ -798,39 +789,20 @@ def decompose(hf_file, mapping="jordan_wigner", core=None, active=None):
 
     mapping = mapping.strip().lower()
 
-    if mapping not in ("jordan_wigner", "bravyi_kitaev"):
-        raise TypeError(
-            f"The '{mapping}' transformation is not available. \n "
-            f"Please set 'mapping' to 'jordan_wigner' or 'bravyi_kitaev'."
-        )
-
     # fermionic-to-qubit transformation of the Hamiltonian
+    if mapping == "parity":
+        binary_code = openfermion.parity_code(molecule.n_qubits)
+        return openfermion.transforms.binary_code_transform(fermionic_hamiltonian, binary_code)
     if mapping == "bravyi_kitaev":
         return openfermion.transforms.bravyi_kitaev(fermionic_hamiltonian)
 
     return openfermion.transforms.jordan_wigner(fermionic_hamiltonian)
 
 
-def molecular_hamiltonian(
-    symbols,
-    coordinates,
-    name="molecule",
-    charge=0,
-    mult=1,
-    basis="sto-3g",
-    method="dhf",
-    active_electrons=None,
-    active_orbitals=None,
-    mapping="jordan_wigner",
-    outpath=".",
-    wires=None,
-    alpha=None,
-    coeff=None,
-    args=None,
-    load_data=False,
-    convert_tol=1e012,
-):  # pylint:disable=too-many-arguments, too-many-statements
-    r"""Generate the qubit Hamiltonian of a molecule.
+def molecular_hamiltonian(*args, **kwargs):
+    """molecular_hamiltonian(molecule, method="dhf", active_electrons=None, active_orbitals=None,\
+    mapping="jordan_wigner", outpath=".", wires=None, args=None, convert_tol=1e12)
+    Generate the qubit Hamiltonian of a molecule.
 
     This function drives the construction of the second-quantized electronic Hamiltonian
     of a molecule and its transformation to the basis of Pauli matrices.
@@ -854,17 +826,7 @@ def molecular_hamiltonian(
     |
 
     Args:
-        symbols (list[str]): symbols of the atomic species in the molecule
-        coordinates (array[float]): atomic positions in Cartesian coordinates.
-            The atomic coordinates must be in atomic units and can be given as either a 1D array of
-            size ``3*N``, or a 2D array of shape ``(N, 3)`` where ``N`` is the number of atoms.
-        name (str): name of the molecule
-        charge (int): Net charge of the molecule. If not specified a neutral system is assumed.
-        mult (int): Spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1`
-            for :math:`N_\mathrm{unpaired}` unpaired electrons occupying the HF orbitals.
-            Possible values of ``mult`` are :math:`1, 2, 3, \ldots`. If not specified,
-            a closed-shell HF state is assumed.
-        basis (str): atomic basis set used to represent the molecular orbitals
+        molecule (~qchem.molecule.Molecule): the molecule object
         method (str): Quantum chemistry method used to solve the
             mean field electronic structure problem. Available options are ``method="dhf"``
             to specify the built-in differentiable Hartree-Fock solver, ``method="pyscf"`` to use
@@ -881,22 +843,26 @@ def molecular_hamiltonian(
             corresponding to the qubit number equal to its index.
             For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted for
             partial mapping. If None, will use identity map.
-        alpha (array[float]): exponents of the primitive Gaussian functions
-        coeff (array[float]): coefficients of the contracted Gaussian functions
         args (array[array[float]]): initial values of the differentiable parameters
-        load_data (bool): flag to load data from the basis-set-exchange library
         convert_tol (float): Tolerance in `machine epsilon <https://numpy.org/doc/stable/reference/generated/numpy.real_if_close.html>`_
             for the imaginary part of the Hamiltonian coefficients created by openfermion.
             Coefficients with imaginary part less than 2.22e-16*tol are considered to be real.
 
+
     Returns:
-        tuple[pennylane.Hamiltonian, int]: the fermionic-to-qubit transformed Hamiltonian
+        tuple[pennylane.Hamiltonian, int]: the fermionic-to-qubit transformed  Hamiltonian
         and the number of qubits
 
+    .. note::
+        The ``molecular_hamiltonian`` function accepts a ``Molecule`` object as its first argument.
+        Look at the `Usage Details` for more details on the old interface.
+
     **Example**
 
-    >>> symbols, coordinates = (['H', 'H'], np.array([0., 0., -0.66140414, 0., 0., 0.66140414]))
-    >>> H, qubits = molecular_hamiltonian(symbols, coordinates)
+    >>> symbols = ['H', 'H']
+    >>> coordinates = np.array([[0., 0., -0.66140414], [0., 0., 0.66140414]])
+    >>> molecule = qml.qchem.Molecule(symbols, coordinates)
+    >>> H, qubits = qml.qchem.molecular_hamiltonian(molecule)
     >>> print(qubits)
     4
     >>> print(H)
@@ -915,11 +881,183 @@ def molecular_hamiltonian(
     + (0.1676831945771896) [Z1 Z2]
     + (0.12293305056183801) [Z1 Z3]
     + (0.176276408043196) [Z2 Z3]
+
+    .. details::
+        :title: Usage Details
+
+        The old interface for this method involved passing molecular information as separate arguments:
+
+            ``molecular_hamiltonian``\\ (`symbols, coordinates, name='molecule', charge=0, mult=1, basis='sto-3g',`
+            `method='dhf', active_electrons=None, active_orbitals=None, mapping='jordan_wigner', outpath='.',`
+            `wires=None, alpha=None, coeff=None, args=None, load_data=False, convert_tol=1e12`)
+
+        Molecule-based Arguments:
+          - **symbols** (list[str]): symbols of the atomic species in the molecule
+          - **coordinates** (array[float]): atomic positions in Cartesian coordinates.
+            The atomic coordinates must be in atomic units and can be given as either a 1D array of
+            size ``3*N``, or a 2D array of shape ``(N, 3)`` where ``N`` is the number of atoms.
+            name (str): name of the molecule
+          - **charge** (int): Net charge of the molecule. If not specified a neutral system is assumed.
+          - **mult** (int): Spin multiplicity :math:`\\mathrm{mult}=N_\\mathrm{unpaired} + 1` for :math:`N_\\mathrm{unpaired}`
+            unpaired electrons occupying the HF orbitals. Possible values of ``mult`` are :math:`1, 2, 3, \\ldots`.
+            If not specified, a closed-shell HF state is assumed.
+          - **basis** (str): atomic basis set used to represent the molecular orbitals
+          - **alpha** (array[float]): exponents of the primitive Gaussian functions
+          - **coeff** (array[float]): coefficients of the contracted Gaussian functions
+
+        Therefore, a molecular Hamiltonian had to be constructed in the following manner:
+
+        .. code-block:: python
+
+            from pennylane import qchem
+
+            symbols = ["H", "H"]
+            geometry = [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
+
+            H, qubit = qchem.molecular_hamiltonian(symbols, geometry, charge=0)
+
+        As part of the new interface, we are shifting towards extracting all the molecular information
+        from the :class:`~.qchem.molecule.Molecule` within the ``molecular_hamiltonian`` method.
+
     """
 
+    if len(args) != 0:
+        return _molecular_hamiltonian_dispatch(*args, **kwargs)
+
+    method = kwargs.pop("symbols", None) or kwargs.pop("molecule", None)
+    if method is not None:
+        return _molecular_hamiltonian_dispatch(method, **kwargs)
+
+    raise NotImplementedError(
+        "The provided arguments do not contain information about symbols in the molecule. "
+        "Please provide that information in the form of a molecule object or as a list of symbols."
+    )
+
+
+@singledispatch
+def _molecular_hamiltonian_dispatch(*args, **kwargs):
+    r"""Generate the qubit Hamiltonian of a molecule."""
+    raise NotImplementedError(
+        "molecular_hamiltonian supports only list or molecule object types. "
+        "Please provide one of them."
+    )
+
+
+@_molecular_hamiltonian_dispatch.register(Molecule)
+def _(
+    molecule,
+    method="dhf",
+    active_electrons=None,
+    active_orbitals=None,
+    mapping="jordan_wigner",
+    outpath=".",
+    wires=None,
+    args=None,
+    convert_tol=1e12,
+):
+    return _molecular_hamiltonian(
+        molecule.symbols,
+        molecule.coordinates,
+        molecule.name,
+        molecule.charge,
+        molecule.mult,
+        molecule.basis_name,
+        method,
+        active_electrons,
+        active_orbitals,
+        mapping,
+        outpath,
+        wires,
+        molecule.alpha,
+        molecule.coeff,
+        args,
+        molecule.load_data,
+        convert_tol,
+    )
+
+
+@_molecular_hamiltonian_dispatch.register(list)
+def _(
+    symbols,
+    coordinates,
+    unit="bohr",
+    name="molecule",
+    charge=0,
+    mult=1,
+    basis="sto-3g",
+    method="dhf",
+    active_electrons=None,
+    active_orbitals=None,
+    mapping="jordan_wigner",
+    outpath=".",
+    wires=None,
+    alpha=None,
+    coeff=None,
+    args=None,
+    load_data=False,
+    convert_tol=1e12,
+):
+
+    if (coord_unit := unit.strip().lower()) not in ("angstrom", "bohr"):
+        raise ValueError(
+            f"The provided unit '{unit}' is not supported. "
+            f"Please set 'unit' to 'bohr' or 'angstrom'."
+        )
+
+    if coord_unit == "angstrom":
+        coordinates = coordinates / bohr_angs
+
+    return _molecular_hamiltonian(
+        symbols,
+        coordinates=coordinates,
+        name=name,
+        charge=charge,
+        mult=mult,
+        basis=basis,
+        method=method,
+        active_electrons=active_electrons,
+        active_orbitals=active_orbitals,
+        mapping=mapping,
+        outpath=outpath,
+        wires=wires,
+        alpha=alpha,
+        coeff=coeff,
+        args=args,
+        load_data=load_data,
+        convert_tol=convert_tol,
+    )
+
+
+def _molecular_hamiltonian(
+    symbols=None,
+    coordinates=None,
+    name="molecule",
+    charge=0,
+    mult=1,
+    basis="sto-3g",
+    method="dhf",
+    active_electrons=None,
+    active_orbitals=None,
+    mapping="jordan_wigner",
+    outpath=".",
+    wires=None,
+    alpha=None,
+    coeff=None,
+    args=None,
+    load_data=False,
+    convert_tol=1e12,
+):  # pylint:disable=too-many-arguments, too-many-statements
+    r"""Generate the qubit Hamiltonian of a molecule."""
+
     if method not in ["dhf", "pyscf", "openfermion"]:
         raise ValueError("Only 'dhf', 'pyscf' and 'openfermion' backends are supported.")
 
+    if mapping.strip().lower() not in ["jordan_wigner", "parity", "bravyi_kitaev"]:
+        raise ValueError(
+            f"'{mapping}' is not supported."
+            f"Please set the mapping to 'jordan_wigner', 'parity' or 'bravyi_kitaev'."
+        )
+
     if len(coordinates) == len(symbols) * 3:
         geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3))
         geometry_hf = coordinates
@@ -933,17 +1071,17 @@ def molecular_hamiltonian(
         wires_new = qml.qchem.convert._process_wires(wires)
         wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels)))
 
-    if method == "dhf":
+    if method in ("dhf", "pyscf"):
+        n_electrons = sum([atomic_numbers[s] for s in symbols]) - charge
 
-        if mapping != "jordan_wigner":
-            raise ValueError(
-                "Only 'jordan_wigner' mapping is supported for the differentiable workflow."
-            )
-        if mult != 1:
+        if n_electrons % 2 == 1 or mult != 1:
             raise ValueError(
-                "Openshell systems are not supported for the differentiable workflow. Use "
-                "`method = 'pyscf'` or change the charge or spin multiplicity of the molecule."
+                "Open-shell systems are not supported for the requested backend. Use "
+                "method = 'openfermion' or change the charge or spin multiplicity of the molecule."
             )
+
+    if method == "dhf":
+
         if args is None and isinstance(geometry_dhf, qml.numpy.tensor):
             geometry_dhf.requires_grad = False
         mol = qml.qchem.Molecule(
@@ -962,9 +1100,9 @@ def molecular_hamiltonian(
 
         requires_grad = args is not None
         h = (
-            qml.qchem.diff_hamiltonian(mol, core=core, active=active)(*args)
+            qml.qchem.diff_hamiltonian(mol, core=core, active=active, mapping=mapping)(*args)
             if requires_grad
-            else qml.qchem.diff_hamiltonian(mol, core=core, active=active)()
+            else qml.qchem.diff_hamiltonian(mol, core=core, active=active, mapping=mapping)()
         )
 
         if active_new_opmath():
@@ -985,19 +1123,33 @@ def molecular_hamiltonian(
             h = qml.map_wires(h, wires_map)
         return h, 2 * len(active)
 
-    if method == "pyscf" and mapping.strip().lower() == "jordan_wigner":
+    if method == "pyscf":
         core_constant, one_mo, two_mo = _pyscf_integrals(
             symbols, geometry_hf, charge, mult, basis, active_electrons, active_orbitals
         )
 
         hf = qml.qchem.fermionic_observable(core_constant, one_mo, two_mo)
+        mapping = mapping.strip().lower()
+        qubits = len(hf.wires)
 
         if active_new_opmath():
-            h_pl = qml.jordan_wigner(hf, wire_map=wires_map, tol=1.0e-10).simplify()
-
+            if mapping == "jordan_wigner":
+                h_pl = qml.jordan_wigner(hf, wire_map=wires_map, tol=1.0e-10)
+            elif mapping == "parity":
+                h_pl = qml.parity_transform(hf, qubits, wire_map=wires_map, tol=1.0e-10)
+            elif mapping == "bravyi_kitaev":
+                h_pl = qml.bravyi_kitaev(hf, qubits, wire_map=wires_map, tol=1.0e-10)
+
+            h_pl.simplify()
         else:
-            h_pl = qml.jordan_wigner(hf, ps=True, wire_map=wires_map, tol=1.0e-10).hamiltonian()
-            h_pl = simplify(h_pl)
+            if mapping == "jordan_wigner":
+                h_pl = qml.jordan_wigner(hf, ps=True, wire_map=wires_map, tol=1.0e-10)
+            elif mapping == "parity":
+                h_pl = qml.parity_transform(hf, qubits, ps=True, wire_map=wires_map, tol=1.0e-10)
+            elif mapping == "bravyi_kitaev":
+                h_pl = qml.bravyi_kitaev(hf, qubits, ps=True, wire_map=wires_map, tol=1.0e-10)
+
+            h_pl = simplify(h_pl.hamiltonian())
 
         return h_pl, len(h_pl.wires)
 
diff --git a/requirements-ci.txt b/requirements-ci.txt
index 6ff2c5afc..932820a26 100644
--- a/requirements-ci.txt
+++ b/requirements-ci.txt
@@ -1,5 +1,5 @@
 numpy
-scipy
+scipy<1.13.0
 cvxpy
 cvxopt
 networkx
@@ -8,7 +8,7 @@ autograd
 toml
 appdirs
 semantic_version
-autoray>=0.6.11
+autoray>=0.6.1,<0.6.10
 matplotlib
 requests
 tomli # Drop once minimum Python version is 3.11


Test patch:
diff --git a/tests/qchem/of_tests/test_decompose.py b/tests/qchem/of_tests/test_decompose.py
index ccb4f4314..88adde1f0 100644
--- a/tests/qchem/of_tests/test_decompose.py
+++ b/tests/qchem/of_tests/test_decompose.py
@@ -539,17 +539,3 @@ def test_transformation(name, core, active, mapping, coeffs_ref, pauli_strings_r
 
     assert np.allclose(coeffs, coeffs_ref, **tol)
     assert pauli_strings == pauli_strings_ref
-
-
-@pytest.mark.usefixtures("skip_if_no_openfermion_support")
-def test_not_available_transformation():
-    r"""Test that an error is raised if the chosen fermionic-to-qubit transformation
-    is neither 'jordan_wigner' nor 'bravyi_kitaev'."""
-
-    with pytest.raises(TypeError, match="transformation is not available"):
-        qchem.decompose(
-            os.path.join(ref_dir, "lih"),
-            mapping="not_available_transformation",
-            core=[0],
-            active=[1, 2],
-        )
diff --git a/tests/qchem/of_tests/test_molecular_hamiltonian.py b/tests/qchem/of_tests/test_molecular_hamiltonian.py
index fc449b3a2..f412f0814 100644
--- a/tests/qchem/of_tests/test_molecular_hamiltonian.py
+++ b/tests/qchem/of_tests/test_molecular_hamiltonian.py
@@ -18,7 +18,7 @@ Unit tests for molecular Hamiltonians.
 import pytest
 
 import pennylane as qml
-from pennylane import Identity, PauliX, PauliY, PauliZ
+from pennylane import I, X, Y, Z
 from pennylane import numpy as np
 from pennylane import qchem
 from pennylane.operation import active_new_opmath
@@ -65,8 +65,8 @@ test_coordinates = np.array(
     ),
     [
         (0, 1, "pyscf", 2, 2, "jordan_WIGNER"),
-        (1, 2, "pyscf", 3, 4, "BRAVYI_kitaev"),
-        (-1, 2, "pyscf", 1, 2, "jordan_WIGNER"),
+        (1, 2, "openfermion", 3, 4, "BRAVYI_kitaev"),
+        (-1, 2, "openfermion", 1, 2, "jordan_WIGNER"),
         (2, 1, "pyscf", 2, 2, "BRAVYI_kitaev"),
     ],
 )
@@ -106,16 +106,65 @@ def test_building_hamiltonian(
 
 
 @pytest.mark.parametrize(
-    ("symbols", "geometry", "h_ref_data"),
+    (
+        "charge",
+        "mult",
+        "package",
+        "nact_els",
+        "nact_orbs",
+        "mapping",
+    ),
+    [
+        (0, 1, "pyscf", 2, 2, "jordan_WIGNER"),
+        (1, 2, "openfermion", 3, 4, "BRAVYI_kitaev"),
+        (-1, 2, "openfermion", 1, 2, "jordan_WIGNER"),
+        (2, 1, "pyscf", 2, 2, "BRAVYI_kitaev"),
+    ],
+)
+@pytest.mark.usefixtures("skip_if_no_openfermion_support", "use_legacy_and_new_opmath")
+def test_building_hamiltonian_molecule_class(
+    charge,
+    mult,
+    package,
+    nact_els,
+    nact_orbs,
+    mapping,
+    tmpdir,
+):
+    r"""Test that the generated Hamiltonian `built_hamiltonian` using the molecule class, is an
+    instance of the PennyLane Hamiltonian class and the correctness of the total number of qubits
+    required to run the quantum simulation. The latter is tested for different values of the
+    molecule's charge and for active spaces with different size"""
+
+    args = qchem.Molecule(test_symbols, test_coordinates, charge=charge, mult=mult)
+    kwargs = {
+        "method": package,
+        "active_electrons": nact_els,
+        "active_orbitals": nact_orbs,
+        "mapping": mapping,
+        "outpath": tmpdir.strpath,
+    }
+
+    built_hamiltonian, qubits = qchem.molecular_hamiltonian(args, **kwargs)
+
+    if active_new_opmath():
+        assert not isinstance(built_hamiltonian, qml.Hamiltonian)
+    else:
+        assert isinstance(built_hamiltonian, qml.Hamiltonian)
+    assert qubits == 2 * nact_orbs
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "mapping", "h_ref_data"),
     [
         (
             ["H", "H"],
             np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+            "jordan_wigner",
             # computed with OpenFermion; data reordered
             # h_mol = molecule.get_molecular_hamiltonian()
             # h_f = openfermion.transforms.get_fermion_operator(h_mol)
             # h_q = openfermion.transforms.jordan_wigner(h_f)
-            # h_pl = qchem.convert_observable(h_q, wires=[0, 1, 2, 3], tol=(5e-5))
             (
                 np.array(
                     [
@@ -137,32 +186,127 @@ def test_building_hamiltonian(
                     ]
                 ),
                 [
-                    Identity(wires=[0]),
-                    PauliZ(wires=[0]),
-                    PauliZ(wires=[1]),
-                    PauliZ(wires=[0]) @ PauliZ(wires=[1]),
-                    PauliY(wires=[0]) @ PauliX(wires=[1]) @ PauliX(wires=[2]) @ PauliY(wires=[3]),
-                    PauliY(wires=[0]) @ PauliY(wires=[1]) @ PauliX(wires=[2]) @ PauliX(wires=[3]),
-                    PauliX(wires=[0]) @ PauliX(wires=[1]) @ PauliY(wires=[2]) @ PauliY(wires=[3]),
-                    PauliX(wires=[0]) @ PauliY(wires=[1]) @ PauliY(wires=[2]) @ PauliX(wires=[3]),
-                    PauliZ(wires=[2]),
-                    PauliZ(wires=[0]) @ PauliZ(wires=[2]),
-                    PauliZ(wires=[3]),
-                    PauliZ(wires=[0]) @ PauliZ(wires=[3]),
-                    PauliZ(wires=[1]) @ PauliZ(wires=[2]),
-                    PauliZ(wires=[1]) @ PauliZ(wires=[3]),
-                    PauliZ(wires=[2]) @ PauliZ(wires=[3]),
+                    I(0),
+                    Z(0),
+                    Z(1),
+                    Z(0) @ Z(1),
+                    Y(0) @ X(1) @ X(2) @ Y(3),
+                    Y(0) @ Y(1) @ X(2) @ X(3),
+                    X(0) @ X(1) @ Y(2) @ Y(3),
+                    X(0) @ Y(1) @ Y(2) @ X(3),
+                    Z(2),
+                    Z(0) @ Z(2),
+                    Z(3),
+                    Z(0) @ Z(3),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
+                ],
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+            "parity",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # binary_code = openfermion.parity_code(molecule.n_qubits)
+            # h_q = openfermion.transforms.binary_code_transform(h_f, binary_code)
+            (
+                np.array(
+                    [
+                        0.2981787007221673,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.20813364101195764,
+                        0.20813364101195767,
+                        0.13290292584331462,
+                        0.13290292584331462,
+                        0.175463287257566,
+                        0.175463287257566,
+                        0.17860976802544348,
+                        -0.34724871015550757,
+                        0.18470917137696227,
+                        -0.3472487101555076,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ Z(1) @ X(2),
+                    X(0) @ Z(1) @ X(2) @ Z(3),
+                    Y(0) @ Y(2),
+                    Y(0) @ Y(2) @ Z(3),
+                    Z(0),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(1) @ Z(2),
+                    Z(0) @ Z(1) @ Z(2) @ Z(3),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(2) @ Z(3),
+                    Z(1),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
+                ],
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+            "bravyi_kitaev",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.bravyi_kitaev(h_f)
+            (
+                np.array(
+                    [
+                        0.2981787007221673,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.20813364101195764,
+                        0.20813364101195767,
+                        0.175463287257566,
+                        0.175463287257566,
+                        0.13290292584331462,
+                        0.13290292584331462,
+                        0.17860976802544348,
+                        -0.3472487101555076,
+                        0.18470917137696227,
+                        -0.34724871015550757,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ Z(1) @ X(2),
+                    X(0) @ Z(1) @ X(2) @ Z(3),
+                    Y(0) @ Z(1) @ Y(2),
+                    Y(0) @ Z(1) @ Y(2) @ Z(3),
+                    Z(0),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(1) @ Z(2),
+                    Z(0) @ Z(1) @ Z(2) @ Z(3),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(2) @ Z(3),
+                    Z(1),
+                    Z(1) @ Z(2) @ Z(3),
+                    Z(1) @ Z(3),
+                    Z(2),
                 ],
             ),
         ),
         (
             ["H", "H"],
             np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+            "jordan_wigner",
             # computed with OpenFermion; data reordered
             # h_mol = molecule.get_molecular_hamiltonian()
             # h_f = openfermion.transforms.get_fermion_operator(h_mol)
             # h_q = openfermion.transforms.jordan_wigner(h_f)
-            # h_pl = qchem.convert_observable(h_q, wires=[0, 1, 2, 3], tol=(5e-5))
             (
                 np.array(
                     [
@@ -184,37 +328,39 @@ def test_building_hamiltonian(
                     ]
                 ),
                 [
-                    Identity(wires=[0]),
-                    PauliZ(wires=[0]),
-                    PauliZ(wires=[1]),
-                    PauliZ(wires=[0]) @ PauliZ(wires=[1]),
-                    PauliY(wires=[0]) @ PauliX(wires=[1]) @ PauliX(wires=[2]) @ PauliY(wires=[3]),
-                    PauliY(wires=[0]) @ PauliY(wires=[1]) @ PauliX(wires=[2]) @ PauliX(wires=[3]),
-                    PauliX(wires=[0]) @ PauliX(wires=[1]) @ PauliY(wires=[2]) @ PauliY(wires=[3]),
-                    PauliX(wires=[0]) @ PauliY(wires=[1]) @ PauliY(wires=[2]) @ PauliX(wires=[3]),
-                    PauliZ(wires=[2]),
-                    PauliZ(wires=[0]) @ PauliZ(wires=[2]),
-                    PauliZ(wires=[3]),
-                    PauliZ(wires=[0]) @ PauliZ(wires=[3]),
-                    PauliZ(wires=[1]) @ PauliZ(wires=[2]),
-                    PauliZ(wires=[1]) @ PauliZ(wires=[3]),
-                    PauliZ(wires=[2]) @ PauliZ(wires=[3]),
+                    I(0),
+                    Z(0),
+                    Z(1),
+                    Z(0) @ Z(1),
+                    Y(0) @ X(1) @ X(2) @ Y(3),
+                    Y(0) @ Y(1) @ X(2) @ X(3),
+                    X(0) @ X(1) @ Y(2) @ Y(3),
+                    X(0) @ Y(1) @ Y(2) @ X(3),
+                    Z(2),
+                    Z(0) @ Z(2),
+                    Z(3),
+                    Z(0) @ Z(3),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
                 ],
             ),
         ),
     ],
 )
 @pytest.mark.usefixtures("use_legacy_and_new_opmath")
-def test_differentiable_hamiltonian(symbols, geometry, h_ref_data):
+def test_differentiable_hamiltonian(symbols, geometry, mapping, h_ref_data):
     r"""Test that molecular_hamiltonian returns the correct Hamiltonian with the differentiable
     backend."""
 
     geometry.requires_grad = True
     args = [geometry.reshape(2, 3)]
-    h_args = qchem.molecular_hamiltonian(symbols, geometry, method="dhf", args=args)[0]
+    h_args = qchem.molecular_hamiltonian(
+        symbols, geometry, method="dhf", args=args, mapping=mapping
+    )[0]
 
     geometry.requires_grad = False
-    h_noargs = qchem.molecular_hamiltonian(symbols, geometry, method="dhf")[0]
+    h_noargs = qchem.molecular_hamiltonian(symbols, geometry, method="dhf", mapping=mapping)[0]
 
     ops = [
         qml.operation.Tensor(*op) if isinstance(op, qml.ops.Prod) else op
@@ -240,28 +386,254 @@ def test_differentiable_hamiltonian(symbols, geometry, h_ref_data):
     )
 
 
-@pytest.mark.usefixtures("use_legacy_and_new_opmath")
 @pytest.mark.parametrize(
-    ("symbols", "geometry", "method", "wiremap"),
+    ("symbols", "geometry", "mapping", "h_ref_data"),
     [
         (
             ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
-            "pyscf",
-            ["a", "b", "c", "d"],
+            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+            "jordan_wigner",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.jordan_wigner(h_f)
+            (
+                np.array(
+                    [
+                        0.2981788017,
+                        0.2081336485,
+                        0.2081336485,
+                        0.1786097698,
+                        0.042560361,
+                        -0.042560361,
+                        -0.042560361,
+                        0.042560361,
+                        -0.3472487379,
+                        0.1329029281,
+                        -0.3472487379,
+                        0.175463289,
+                        0.175463289,
+                        0.1329029281,
+                        0.1847091733,
+                    ]
+                ),
+                [
+                    I(0),
+                    Z(0),
+                    Z(1),
+                    Z(0) @ Z(1),
+                    Y(0) @ X(1) @ X(2) @ Y(3),
+                    Y(0) @ Y(1) @ X(2) @ X(3),
+                    X(0) @ X(1) @ Y(2) @ Y(3),
+                    X(0) @ Y(1) @ Y(2) @ X(3),
+                    Z(2),
+                    Z(0) @ Z(2),
+                    Z(3),
+                    Z(0) @ Z(3),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
+                ],
+            ),
         ),
         (
             ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
-            "pyscf",
-            [0, "z", 3, "ancilla"],
+            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+            "parity",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # binary_code = openfermion.parity_code(molecule.n_qubits)
+            # h_q = openfermion.transforms.binary_code_transform(h_f, binary_code)
+            (
+                np.array(
+                    [
+                        0.2981787007221673,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.20813364101195764,
+                        0.20813364101195767,
+                        0.13290292584331462,
+                        0.13290292584331462,
+                        0.175463287257566,
+                        0.175463287257566,
+                        0.17860976802544348,
+                        -0.34724871015550757,
+                        0.18470917137696227,
+                        -0.3472487101555076,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ Z(1) @ X(2),
+                    X(0) @ Z(1) @ X(2) @ Z(3),
+                    Y(0) @ Y(2),
+                    Y(0) @ Y(2) @ Z(3),
+                    Z(0),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(1) @ Z(2),
+                    Z(0) @ Z(1) @ Z(2) @ Z(3),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(2) @ Z(3),
+                    Z(1),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
+                ],
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+            "bravyi_kitaev",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.bravyi_kitaev(h_f)
+            (
+                np.array(
+                    [
+                        0.2981787007221673,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.04256036141425139,
+                        0.20813364101195764,
+                        0.20813364101195767,
+                        0.175463287257566,
+                        0.175463287257566,
+                        0.13290292584331462,
+                        0.13290292584331462,
+                        0.17860976802544348,
+                        -0.3472487101555076,
+                        0.18470917137696227,
+                        -0.34724871015550757,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ Z(1) @ X(2),
+                    X(0) @ Z(1) @ X(2) @ Z(3),
+                    Y(0) @ Z(1) @ Y(2),
+                    Y(0) @ Z(1) @ Y(2) @ Z(3),
+                    Z(0),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(1) @ Z(2),
+                    Z(0) @ Z(1) @ Z(2) @ Z(3),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(2) @ Z(3),
+                    Z(1),
+                    Z(1) @ Z(2) @ Z(3),
+                    Z(1) @ Z(3),
+                    Z(2),
+                ],
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]),
+            "jordan_wigner",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.jordan_wigner(h_f)
+            (
+                np.array(
+                    [
+                        0.2981788017,
+                        0.2081336485,
+                        0.2081336485,
+                        0.1786097698,
+                        0.042560361,
+                        -0.042560361,
+                        -0.042560361,
+                        0.042560361,
+                        -0.3472487379,
+                        0.1329029281,
+                        -0.3472487379,
+                        0.175463289,
+                        0.175463289,
+                        0.1329029281,
+                        0.1847091733,
+                    ]
+                ),
+                [
+                    I(0),
+                    Z(0),
+                    Z(1),
+                    Z(0) @ Z(1),
+                    Y(0) @ X(1) @ X(2) @ Y(3),
+                    Y(0) @ Y(1) @ X(2) @ X(3),
+                    X(0) @ X(1) @ Y(2) @ Y(3),
+                    X(0) @ Y(1) @ Y(2) @ X(3),
+                    Z(2),
+                    Z(0) @ Z(2),
+                    Z(3),
+                    Z(0) @ Z(3),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
+                ],
+            ),
         ),
     ],
 )
+@pytest.mark.usefixtures("use_legacy_and_new_opmath")
+def test_differentiable_hamiltonian_molecule_class(symbols, geometry, mapping, h_ref_data):
+    r"""Test that molecular_hamiltonian generated using the molecule class
+    returns the correct Hamiltonian with the differentiable backend."""
+
+    geometry.requires_grad = True
+    args = [geometry.reshape(2, 3)]
+    molecule = qchem.Molecule(symbols, geometry)
+    h_args = qchem.molecular_hamiltonian(molecule, method="dhf", args=args, mapping=mapping)[0]
+
+    geometry.requires_grad = False
+    molecule = qchem.Molecule(symbols, geometry)
+    h_noargs = qchem.molecular_hamiltonian(molecule, method="dhf", mapping=mapping)[0]
+
+    ops = [
+        qml.operation.Tensor(*op) if isinstance(op, qml.ops.Prod) else op
+        for op in map(qml.simplify, h_ref_data[1])
+    ]
+    h_ref = qml.Hamiltonian(h_ref_data[0], ops)
+
+    h_ref_coeffs, h_ref_ops = h_ref.terms()
+    h_args_coeffs, h_args_ops = h_args.terms()
+    h_noargs_coeffs, h_noargs_ops = h_noargs.terms()
+
+    assert all(coeff.requires_grad is True for coeff in h_args_coeffs)
+    assert all(coeff.requires_grad is False for coeff in h_noargs_coeffs)
+
+    assert np.allclose(np.sort(h_args_coeffs), np.sort(h_ref_coeffs))
+    assert qml.Hamiltonian(np.ones(len(h_args_coeffs)), h_args_ops).compare(
+        qml.Hamiltonian(np.ones(len(h_ref_coeffs)), h_ref_ops)
+    )
+
+    assert np.allclose(np.sort(h_noargs_coeffs), np.sort(h_ref_coeffs))
+    assert qml.Hamiltonian(np.ones(len(h_noargs_coeffs)), h_noargs_ops).compare(
+        qml.Hamiltonian(np.ones(len(h_ref_coeffs)), h_ref_ops)
+    )
+
+
+@pytest.mark.usefixtures("use_legacy_and_new_opmath")
+@pytest.mark.parametrize(
+    ("wiremap"),
+    [
+        ["a", "b", "c", "d"],
+        [0, "z", 3, "ancilla"],
+    ],
+)
 @pytest.mark.usefixtures("skip_if_no_openfermion_support")
-def test_custom_wiremap_hamiltonian_pyscf(symbols, geometry, method, wiremap, tmpdir):
+def test_custom_wiremap_hamiltonian_pyscf(wiremap, tmpdir):
     r"""Test that the generated Hamiltonian has the correct wire labels given by a custom wiremap."""
 
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
+    method = "pyscf"
+
     hamiltonian, _ = qchem.molecular_hamiltonian(
         symbols=symbols,
         coordinates=geometry,
@@ -273,27 +645,50 @@ def test_custom_wiremap_hamiltonian_pyscf(symbols, geometry, method, wiremap, tm
     assert set(hamiltonian.wires) == set(wiremap)
 
 
+@pytest.mark.parametrize(
+    ("wiremap"),
+    [
+        ["a", "b", "c", "d"],
+        [0, "z", 3, "ancilla"],
+    ],
+)
+@pytest.mark.usefixtures("skip_if_no_openfermion_support")
+def test_custom_wiremap_hamiltonian_pyscf_molecule_class(wiremap, tmpdir):
+    r"""Test that the generated Hamiltonian has the correct wire labels given by a custom wiremap."""
+
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
+    method = "pyscf"
+    molecule = qchem.Molecule(symbols, geometry)
+    hamiltonian, _ = qchem.molecular_hamiltonian(
+        molecule,
+        method=method,
+        wires=wiremap,
+        outpath=tmpdir.strpath,
+    )
+
+    assert set(hamiltonian.wires) == set(wiremap)
+
+
 @pytest.mark.usefixtures("use_legacy_and_new_opmath")
 @pytest.mark.parametrize(
-    ("symbols", "geometry", "wiremap", "args"),
+    ("wiremap", "args"),
     [
         (
-            ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
             [0, "z", 3, "ancilla"],
             None,
         ),
         (
-            ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
             [0, "z", 3, "ancilla"],
             [np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])],
         ),
     ],
 )
-def test_custom_wiremap_hamiltonian_dhf(symbols, geometry, wiremap, args, tmpdir):
+def test_custom_wiremap_hamiltonian_dhf(wiremap, args, tmpdir):
     r"""Test that the generated Hamiltonian has the correct wire labels given by a custom wiremap."""
 
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
     wiremap_dict = dict(zip(range(len(wiremap)), wiremap))
 
     hamiltonian_ref, _ = qchem.molecular_hamiltonian(
@@ -316,6 +711,46 @@ def test_custom_wiremap_hamiltonian_dhf(symbols, geometry, wiremap, args, tmpdir
     assert wiremap_calc == wiremap_dict
 
 
+@pytest.mark.usefixtures("use_legacy_and_new_opmath")
+@pytest.mark.parametrize(
+    ("wiremap", "args"),
+    [
+        (
+            [0, "z", 3, "ancilla"],
+            None,
+        ),
+        (
+            [0, "z", 3, "ancilla"],
+            [np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])],
+        ),
+    ],
+)
+def test_custom_wiremap_hamiltonian_dhf_molecule_class(wiremap, args, tmpdir):
+    r"""Test that the generated Hamiltonian has the correct wire labels given by a custom wiremap."""
+
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
+    wiremap_dict = dict(zip(range(len(wiremap)), wiremap))
+
+    molecule = qchem.Molecule(symbols, geometry)
+    hamiltonian_ref, _ = qchem.molecular_hamiltonian(
+        molecule,
+        args=args,
+        outpath=tmpdir.strpath,
+    )
+
+    hamiltonian, _ = qchem.molecular_hamiltonian(
+        molecule,
+        wires=wiremap,
+        args=args,
+        outpath=tmpdir.strpath,
+    )
+
+    wiremap_calc = dict(zip(list(hamiltonian_ref.wires), list(hamiltonian.wires)))
+
+    assert wiremap_calc == wiremap_dict
+
+
 file_content = """\
 2
 in Angstrom
@@ -339,57 +774,92 @@ def test_mol_hamiltonian_with_read_structure(tmpdir):
     assert num_qubits == 4
 
 
-@pytest.mark.parametrize(
-    ("symbols", "geometry"),
-    [
-        (
-            ["H", "H"],
-            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
-        ),
-    ],
-)
-def test_diff_hamiltonian_error(symbols, geometry):
+def test_mol_hamiltonian_with_read_structure_molecule_class(tmpdir):
+    """Test that the pipeline of using molecular_hamiltonian with
+    read_structure executes without errors."""
+    f_name = "h2.xyz"
+    filename = tmpdir.join(f_name)
+
+    with open(filename, "w") as f:
+        f.write(file_content)
+
+    symbols, coordinates = qchem.read_structure(str(filename), outpath=tmpdir)
+
+    molecule = qchem.Molecule(symbols, coordinates)
+    H, num_qubits = qchem.molecular_hamiltonian(molecule)
+    assert len(H.terms()) == 2
+    assert num_qubits == 4
+
+
+def test_diff_hamiltonian_error():
     r"""Test that molecular_hamiltonian raises an error with unsupported mapping."""
 
-    with pytest.raises(ValueError, match="Only 'jordan_wigner' mapping is supported"):
-        qchem.molecular_hamiltonian(symbols, geometry, method="dhf", mapping="bravyi_kitaev")
+    symbols = ["H", "H"]
+    geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0])
 
     with pytest.raises(
         ValueError, match="Only 'dhf', 'pyscf' and 'openfermion' backends are supported"
     ):
         qchem.molecular_hamiltonian(symbols, geometry, method="psi4")
 
-    with pytest.raises(ValueError, match="Openshell systems are not supported"):
+    with pytest.raises(ValueError, match="Open-shell systems are not supported"):
         qchem.molecular_hamiltonian(symbols, geometry, mult=3)
 
 
+def test_pyscf_hamiltonian_error():
+    r"""Test that molecular_hamiltonian raises an error for open-shell systems."""
+
+    symbols = ["H", "H"]
+    geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0])
+
+    with pytest.raises(ValueError, match="Open-shell systems are not supported"):
+        qchem.molecular_hamiltonian(symbols, geometry, mult=3, method="pyscf")
+
+    molecule = qchem.Molecule(symbols, geometry, mult=3)
+    with pytest.raises(ValueError, match="Open-shell systems are not supported"):
+        qchem.molecular_hamiltonian(molecule, method="pyscf")
+
+
+def test_diff_hamiltonian_error_molecule_class():
+    r"""Test that molecular_hamiltonian raises an error with unsupported mapping."""
+
+    symbols = ["H", "H"]
+    geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0])
+
+    molecule = qchem.Molecule(symbols, geometry)
+    with pytest.raises(
+        ValueError, match="Only 'dhf', 'pyscf' and 'openfermion' backends are supported"
+    ):
+        qchem.molecular_hamiltonian(molecule, method="psi4")
+
+    with pytest.raises(ValueError, match="'bksf' is not supported."):
+        qchem.molecular_hamiltonian(molecule, mapping="bksf")
+
+
 @pytest.mark.parametrize(
-    ("symbols", "geometry", "method", "args"),
+    ("method", "args"),
     [
         (
-            ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
             "pyscf",
             None,
         ),
         (
-            ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
             "dhf",
             None,
         ),
         (
-            ["H", "H"],
-            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
             "dhf",
             [np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])],
         ),
     ],
 )
 @pytest.mark.usefixtures("skip_if_no_openfermion_support", "use_legacy_and_new_opmath")
-def test_real_hamiltonian(symbols, geometry, method, args, tmpdir):
+def test_real_hamiltonian(method, args, tmpdir):
     r"""Test that the generated Hamiltonian has real coefficients."""
 
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
+
     hamiltonian, _ = qchem.molecular_hamiltonian(
         symbols=symbols,
         coordinates=geometry,
@@ -401,6 +871,41 @@ def test_real_hamiltonian(symbols, geometry, method, args, tmpdir):
     assert np.isrealobj(hamiltonian.terms()[0])
 
 
+@pytest.mark.parametrize(
+    ("method", "args"),
+    [
+        (
+            "pyscf",
+            None,
+        ),
+        (
+            "dhf",
+            None,
+        ),
+        (
+            "dhf",
+            [np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])],
+        ),
+    ],
+)
+@pytest.mark.usefixtures("skip_if_no_openfermion_support", "use_legacy_and_new_opmath")
+def test_real_hamiltonian_molecule_class(method, args, tmpdir):
+    r"""Test that the generated Hamiltonian has real coefficients."""
+
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
+
+    molecule = qchem.Molecule(symbols, geometry)
+    hamiltonian, _ = qchem.molecular_hamiltonian(
+        molecule,
+        method=method,
+        args=args,
+        outpath=tmpdir.strpath,
+    )
+
+    assert np.isrealobj(hamiltonian.terms()[0])
+
+
 @pytest.mark.parametrize(
     ("symbols", "geometry", "core_ref", "one_ref", "two_ref"),
     [
@@ -433,3 +938,399 @@ def test_pyscf_integrals(symbols, geometry, core_ref, one_ref, two_ref):
     assert np.allclose(core, core_ref)
     assert np.allclose(one, one_ref)
     assert np.allclose(two, two_ref)
+
+
+@pytest.mark.usefixtures("skip_if_no_openfermion_support", "use_legacy_and_new_opmath")
+def test_molecule_as_kwargs(tmpdir):
+    r"""Test that molecular_hamiltonian function works with molecule as
+    keyword argument
+    """
+
+    molecule = qchem.Molecule(
+        test_symbols,
+        test_coordinates,
+    )
+    built_hamiltonian, qubits = qchem.molecular_hamiltonian(
+        molecule=molecule,
+        method="pyscf",
+        active_electrons=2,
+        active_orbitals=2,
+        outpath=tmpdir.strpath,
+    )
+
+    if active_new_opmath():
+        assert not isinstance(built_hamiltonian, qml.Hamiltonian)
+    else:
+        assert isinstance(built_hamiltonian, qml.Hamiltonian)
+    assert qubits == 4
+
+
+def test_error_raised_for_incompatible_type():
+    r"""Test that molecular_hamiltonian raises an error when input is not
+    a list or molecule object.
+    """
+
+    with pytest.raises(
+        NotImplementedError,
+        match="molecular_hamiltonian supports only list or molecule object types.",
+    ):
+        qchem.molecular_hamiltonian(symbols=1, coordinates=test_coordinates, method="dhf")
+
+
+def test_error_raised_for_missing_molecule_information():
+    r"""Test that molecular_hamiltonian raises an error when symbols, and coordinates
+    information is not provided.
+    """
+
+    with pytest.raises(
+        NotImplementedError,
+        match="The provided arguments do not contain information about symbols in the molecule.",
+    ):
+        qchem.molecular_hamiltonian(charge=0, mult=1, method="dhf")
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "charge", "mapping", "h_ref_data"),
+    [
+        (
+            ["H", "H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]]),
+            1,
+            "jordan_wigner",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.jordan_wigner(h_f)
+            (
+                np.array(
+                    [
+                        1.3657458030310135,
+                        -0.03586487568545097,
+                        -0.03201092703651771,
+                        0.03586487568545097,
+                        0.03201092703651771,
+                        -0.031492075818254375,
+                        -0.031492075818254375,
+                        0.037654484403957786,
+                        -0.032229274210852504,
+                        -0.0022222066814484463,
+                        -0.03371428249970282,
+                        -0.030511339285364605,
+                        0.03586487568545097,
+                        0.03201092703651771,
+                        -0.03586487568545097,
+                        -0.03201092703651771,
+                        -0.031492075818254375,
+                        -0.031492075818254375,
+                        0.037654484403957786,
+                        -0.032229274210852504,
+                        -0.0022222066814484463,
+                        -0.03371428249970282,
+                        -0.030511339285364605,
+                        0.27235785388149386,
+                        -0.03051133928536461,
+                        -0.03051133928536461,
+                        0.17448913735995256,
+                        0.11784682872956924,
+                        0.15371170441502022,
+                        0.1487316290904712,
+                        0.18074255612698886,
+                        0.031492075818254375,
+                        -0.031492075818254375,
+                        0.037654484403957786,
+                        -0.032229274210852504,
+                        -0.03371428249970282,
+                        -0.0022222066814484463,
+                        -0.031492075818254375,
+                        0.031492075818254375,
+                        0.037654484403957786,
+                        -0.032229274210852504,
+                        -0.03371428249970282,
+                        -0.0022222066814484463,
+                        0.27235785388149386,
+                        0.15371170441502022,
+                        0.11784682872956924,
+                        0.18074255612698886,
+                        0.1487316290904712,
+                        -0.03583418633226662,
+                        0.03583418633226662,
+                        0.03583418633226662,
+                        -0.03583418633226662,
+                        -0.06458411201474276,
+                        0.16096866344343394,
+                        0.1288375790750158,
+                        0.16467176540728246,
+                        -0.06458411201474279,
+                        0.16467176540728246,
+                        0.1288375790750158,
+                        -0.8044935587718376,
+                        0.20315172438516313,
+                        -0.8044935587718377,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ X(1) @ Y(2) @ Y(3),
+                    X(0) @ X(1) @ Y(4) @ Y(5),
+                    X(0) @ Y(1) @ Y(2) @ X(3),
+                    X(0) @ Y(1) @ Y(4) @ X(5),
+                    X(0) @ Z(1) @ X(2) @ X(3) @ Z(4) @ X(5),
+                    X(0) @ Z(1) @ X(2) @ Y(3) @ Z(4) @ Y(5),
+                    X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4),
+                    X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4) @ Z(5),
+                    X(0) @ Z(1) @ Z(2) @ X(4),
+                    X(0) @ Z(1) @ Z(3) @ X(4),
+                    X(0) @ Z(2) @ Z(3) @ X(4),
+                    Y(0) @ X(1) @ X(2) @ Y(3),
+                    Y(0) @ X(1) @ X(4) @ Y(5),
+                    Y(0) @ Y(1) @ X(2) @ X(3),
+                    Y(0) @ Y(1) @ X(4) @ X(5),
+                    Y(0) @ Z(1) @ Y(2) @ X(3) @ Z(4) @ X(5),
+                    Y(0) @ Z(1) @ Y(2) @ Y(3) @ Z(4) @ Y(5),
+                    Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4),
+                    Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4) @ Z(5),
+                    Y(0) @ Z(1) @ Z(2) @ Y(4),
+                    Y(0) @ Z(1) @ Z(3) @ Y(4),
+                    Y(0) @ Z(2) @ Z(3) @ Y(4),
+                    Z(0),
+                    Z(0) @ X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5),
+                    Z(0) @ Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(3),
+                    Z(0) @ Z(4),
+                    Z(0) @ Z(5),
+                    X(1) @ X(2) @ Y(3) @ Y(4),
+                    X(1) @ Y(2) @ Y(3) @ X(4),
+                    X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5),
+                    X(1) @ Z(2) @ Z(3) @ X(5),
+                    X(1) @ Z(2) @ Z(4) @ X(5),
+                    X(1) @ Z(3) @ Z(4) @ X(5),
+                    Y(1) @ X(2) @ X(3) @ Y(4),
+                    Y(1) @ Y(2) @ X(3) @ X(4),
+                    Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5),
+                    Y(1) @ Z(2) @ Z(3) @ Y(5),
+                    Y(1) @ Z(2) @ Z(4) @ Y(5),
+                    Y(1) @ Z(3) @ Z(4) @ Y(5),
+                    Z(1),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(1) @ Z(4),
+                    Z(1) @ Z(5),
+                    X(2) @ X(3) @ Y(4) @ Y(5),
+                    X(2) @ Y(3) @ Y(4) @ X(5),
+                    Y(2) @ X(3) @ X(4) @ Y(5),
+                    Y(2) @ Y(3) @ X(4) @ X(5),
+                    Z(2),
+                    Z(2) @ Z(3),
+                    Z(2) @ Z(4),
+                    Z(2) @ Z(5),
+                    Z(3),
+                    Z(3) @ Z(4),
+                    Z(3) @ Z(5),
+                    Z(4),
+                    Z(4) @ Z(5),
+                    Z(5),
+                ],
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
+            0,
+            "parity",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # binary_code = openfermion.parity_code(molecule.n_qubits)
+            # h_q = openfermion.transforms.binary_code_transform(h_f, binary_code)
+            (
+                np.array(
+                    [
+                        -0.3596823978788041,
+                        0.050130618654510024,
+                        0.050130618654510024,
+                        0.050130618654510024,
+                        0.050130618654510024,
+                        0.13082413502487947,
+                        0.13082413502487947,
+                        0.1031689785681825,
+                        0.1031689785681825,
+                        0.15329959722269254,
+                        0.15329959722269254,
+                        0.15405495529252655,
+                        -0.11496333923452409,
+                        0.16096866344343408,
+                        -0.11496333923452409,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ Z(1) @ X(2),
+                    X(0) @ Z(1) @ X(2) @ Z(3),
+                    Y(0) @ Y(2),
+                    Y(0) @ Y(2) @ Z(3),
+                    Z(0),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(1) @ Z(2),
+                    Z(0) @ Z(1) @ Z(2) @ Z(3),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(2) @ Z(3),
+                    Z(1),
+                    Z(1) @ Z(2),
+                    Z(1) @ Z(3),
+                    Z(2) @ Z(3),
+                ],
+            ),
+        ),
+        (
+            ["H", "H"],
+            np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]),
+            0,
+            "bravyi_kitaev",
+            # computed with OpenFermion; data reordered
+            # h_mol = molecule.get_molecular_hamiltonian()
+            # h_f = openfermion.transforms.get_fermion_operator(h_mol)
+            # h_q = openfermion.transforms.bravyi_kitaev(h_f)
+            (
+                np.array(
+                    [
+                        -0.3596823978788041,
+                        0.050130618654510024,
+                        0.050130618654510024,
+                        0.050130618654510024,
+                        0.050130618654510024,
+                        0.13082413502487947,
+                        0.13082413502487947,
+                        0.15329959722269254,
+                        0.15329959722269254,
+                        0.1031689785681825,
+                        0.1031689785681825,
+                        0.15405495529252655,
+                        -0.11496333923452409,
+                        0.16096866344343408,
+                        -0.11496333923452409,
+                    ]
+                ),
+                [
+                    I(0),
+                    X(0) @ Z(1) @ X(2),
+                    X(0) @ Z(1) @ X(2) @ Z(3),
+                    Y(0) @ Z(1) @ Y(2),
+                    Y(0) @ Z(1) @ Y(2) @ Z(3),
+                    Z(0),
+                    Z(0) @ Z(1),
+                    Z(0) @ Z(1) @ Z(2),
+                    Z(0) @ Z(1) @ Z(2) @ Z(3),
+                    Z(0) @ Z(2),
+                    Z(0) @ Z(2) @ Z(3),
+                    Z(1),
+                    Z(1) @ Z(2) @ Z(3),
+                    Z(1) @ Z(3),
+                    Z(2),
+                ],
+            ),
+        ),
+    ],
+)
+@pytest.mark.usefixtures("use_legacy_and_new_opmath")
+def test_mapped_hamiltonian_pyscf_openfermion(
+    symbols, geometry, charge, mapping, h_ref_data, tmpdir
+):
+    r"""Test that molecular_hamiltonian returns the correct qubit Hamiltonian with the pyscf and openfermion
+    backend."""
+    methods = ["openfermion", "pyscf"]
+    for method in methods:
+        geometry.requires_grad = False
+        molecule = qchem.Molecule(symbols, geometry, charge=charge)
+        h = qchem.molecular_hamiltonian(
+            molecule, method=method, mapping=mapping, outpath=tmpdir.strpath
+        )[0]
+
+        ops = [
+            qml.operation.Tensor(*op) if isinstance(op, qml.ops.Prod) else op
+            for op in map(qml.simplify, h_ref_data[1])
+        ]
+        h_ref = qml.Hamiltonian(h_ref_data[0], ops)
+
+        h_ref_coeffs, h_ref_ops = h_ref.terms()
+        h_coeffs, h_ops = h.terms()
+
+        assert np.allclose(np.sort(h_coeffs), np.sort(h_ref_coeffs))
+        assert qml.Hamiltonian(np.ones(len(h_coeffs)), h_ops).compare(
+            qml.Hamiltonian(np.ones(len(h_ref_coeffs)), h_ref_ops)
+        )
+
+
+@pytest.mark.parametrize(
+    ("method"),
+    [
+        "pyscf",
+        "dhf",
+        "openfermion",
+    ],
+)
+def test_coordinate_units_for_molecular_hamiltonian(method, tmpdir):
+    r"""Test that molecular_hamiltonian generates the Hamiltonian for both Bohr and Angstrom units."""
+
+    symbols = ["H", "H"]
+    geometry_bohr = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+    geometry_ang = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.529177210903]])
+
+    hamiltonian_bohr, _ = qchem.molecular_hamiltonian(
+        symbols,
+        geometry_bohr,
+        unit="bohr",
+        method=method,
+        outpath=tmpdir.strpath,
+    )
+
+    hamiltonian_ang, _ = qchem.molecular_hamiltonian(
+        symbols,
+        geometry_ang,
+        unit="angstrom",
+        method=method,
+        outpath=tmpdir.strpath,
+    )
+    assert qml.equal(hamiltonian_ang, hamiltonian_bohr)
+
+
+@pytest.mark.parametrize(
+    ("method"),
+    [
+        "pyscf",
+        "dhf",
+        "openfermion",
+    ],
+)
+def test_coordinate_units_for_molecular_hamiltonian_molecule_class(method, tmpdir):
+    r"""Test that molecular_hamiltonian generates the Hamiltonian for both Bohr and Angstrom units."""
+
+    symbols = ["H", "H"]
+    geometry_bohr = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+    geometry_ang = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.529177210903]])
+
+    molecule_bohr = qchem.Molecule(symbols, geometry_bohr, unit="bohr")
+    hamiltonian_bohr, _ = qchem.molecular_hamiltonian(
+        molecule_bohr,
+        method=method,
+        outpath=tmpdir.strpath,
+    )
+
+    molecule_ang = qchem.Molecule(symbols, geometry_ang, unit="angstrom")
+    hamiltonian_ang, _ = qchem.molecular_hamiltonian(
+        molecule_ang,
+        method=method,
+        outpath=tmpdir.strpath,
+    )
+    assert qml.equal(hamiltonian_ang, hamiltonian_bohr)
+
+
+def test_unit_error_molecular_hamiltonian():
+    r"""Test that an error is raised if a wrong/not-supported unit for coordinates is entered."""
+
+    symbols = ["H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+
+    with pytest.raises(ValueError, match="The provided unit 'degrees' is not supported."):
+        qchem.molecular_hamiltonian(symbols, geometry, unit="degrees")
diff --git a/tests/qchem/test_hartree_fock.py b/tests/qchem/test_hartree_fock.py
index a0bcf2b56..2709a528a 100644
--- a/tests/qchem/test_hartree_fock.py
+++ b/tests/qchem/test_hartree_fock.py
@@ -79,6 +79,16 @@ def test_scf(symbols, geometry, v_fock, coeffs, fock_matrix, h_core, repulsion_t
     assert np.allclose(e, repulsion_tensor)
 
 
+def test_scf_openshell_error():
+    r"""Test that scf raises an error when an open-shell molecule is provided."""
+    symbols = ["H", "H", "H"]
+    geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, 2.0]], requires_grad=False)
+    mol = qchem.Molecule(symbols, geometry)
+
+    with pytest.raises(ValueError, match="Open-shell systems are not supported."):
+        qchem.scf(mol)()
+
+
 @pytest.mark.parametrize(
     ("symbols", "geometry", "charge", "basis_name", "e_ref"),
     [
diff --git a/tests/qchem/test_molecule.py b/tests/qchem/test_molecule.py
index 5191f3445..b4e797e56 100644
--- a/tests/qchem/test_molecule.py
+++ b/tests/qchem/test_molecule.py
@@ -46,17 +46,6 @@ class TestMolecule:
         with pytest.raises(ValueError, match="Currently, the only supported basis sets"):
             qchem.Molecule(symbols, geometry, basis_name="6-3_1_g")
 
-    @pytest.mark.parametrize(
-        ("symbols", "geometry", "charge", "mult"),
-        [
-            (["H", "He"], np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]), 0, 2),
-        ],
-    )
-    def test_openshell_error(self, symbols, geometry, charge, mult):
-        r"""Test that an error is raised if the molecule has unpaired electrons."""
-        with pytest.raises(ValueError, match="Openshell systems are not supported"):
-            qchem.Molecule(symbols, geometry, charge=charge, mult=mult)
-
     @pytest.mark.parametrize(
         ("symbols", "geometry"),
         [
@@ -353,3 +342,12 @@ class TestMolecule:
         assert l == l_ref
         assert alpha == alpha_ref
         assert coeff == coeff_ref
+
+    def test_unit_error(self):
+        r"""Test that an error is raised if a wrong/not-supported unit for coordinates is entered."""
+
+        symbols = ["H", "H"]
+        geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+
+        with pytest.raises(ValueError, match="The provided unit 'degrees' is not supported."):
+            qchem.Molecule(symbols, geometry, unit="degrees")