New branch for qchem features (#5638)

**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:**

---------

Co-authored-by: soranjh <soran.jahangiri@gmail.com>
Co-authored-by: Utkarsh <utkarshazad98@gmail.com>
Co-authored-by: Austin Huang <65315367+austingmhuang@users.noreply.github.com>
Co-authored-by: soranjh <40344468+soranjh@users.noreply.github.com>
Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com>

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

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):

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,

pennylane/qchem/hamiltonian.py

@@ -181,14 +181,16 @@ def _fermionic_hamiltonian(*args):
     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:
         mol (~qchem.molecule.Molecule): the molecule object
         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 _molecular_hamiltonian(*args):
 
         h_ferm = fermionic_hamiltonian(mol, cutoff, core, active)(*args)
 
-        return qubit_observable(h_ferm)
+        return qubit_observable(h_ferm, mapping=mapping)
 
     return _molecular_hamiltonian

pennylane/qchem/hartree_fock.py

@@ -116,6 +116,11 @@ def _scf(*args):
             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

pennylane/qchem/molecule.py

@@ -27,6 +27,9 @@
 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 @@ def __init__(
         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 @@ def __init__(
 
         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]
 

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():