Repository: PennyLaneAI/pennylane
Pull request #5118 (2024-01-29 18:54:19)

Add direct call to pyscf in qchem
**Context:**
The `'pyscf'` backend in `qchem.molecular_hamiltonian` currently runs `'OpenFermion-PySCF'` to compute molecular integrals with `PySCF`, and uses these integrals to construct a qubit Hamiltonian with `OpenFermion`. With the new support for fermionic operations in PennyLane, `qchem.molecular_hamiltonian` can directly run `PySCF` when the `'pyscf'` blackened is selected. 

**Description of the Change:**

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**


soranjh: [sc-52875]
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5118?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: `3 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`b5f302c`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/b5f302c5fc903762ea083b5f7df19148cc06c23b?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.69% compared to head [(`3efc09c`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5118?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.67%.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5118?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/qchem/openfermion\_obs.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5118?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3FjaGVtL29wZW5mZXJtaW9uX29icy5weQ==) | 95.08% | [3 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5118?src=pr&el=tree&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    #5118      +/-   ##
==========================================
- Coverage   99.69%   99.67%   -0.02%     
==========================================
  Files         394      394              
  Lines       35979    35770     -209     
==========================================
- Hits        35868    35655     -213     
- Misses        111      115       +4     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5118?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).

trbromley: @soranjh, would you mind adding an entry to the changelog?

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index f4491b31d..ff6e83122 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -6,6 +6,9 @@
 
 <h3>Improvements 🛠</h3>
 
+* The `molecular_hamiltonian` function calls `PySCF` directly when `method='pyscf'` is selected.
+  [(#5118)](https://github.com/PennyLaneAI/pennylane/pull/5118)
+
 <h3>Breaking changes 💔</h3>
 
 <h3>Deprecations 👋</h3>
@@ -21,4 +24,5 @@
 
 This release contains contributions from (in alphabetical order):
 
-Korbinian Kottmann
\ No newline at end of file
+Soran Jahangiri,
+Korbinian Kottmann.
diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py
index 5fa86d490..d7fba9e7f 100644
--- a/pennylane/qchem/openfermion_obs.py
+++ b/pennylane/qchem/openfermion_obs.py
@@ -21,6 +21,8 @@ import numpy as np
 
 import pennylane as qml
 from pennylane.operation import active_new_opmath
+from pennylane.pauli.utils import simplify
+
 
 # Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0)
 bohr_angs = 0.529177210903
@@ -40,6 +42,19 @@ def _import_of():
     return openfermion, openfermionpyscf
 
 
+def _import_pyscf():
+    """Import pyscf."""
+    try:
+        # pylint: disable=import-outside-toplevel, unused-import, multiple-imports
+        import pyscf
+    except ImportError as Error:
+        raise ImportError(
+            "This feature requires pyscf. It can be installed with: pip install pyscf."
+        ) from Error
+
+    return pyscf
+
+
 def observable(fermion_ops, init_term=0, mapping="jordan_wigner", wires=None):
     r"""Builds the fermionic many-body observable whose expectation value can be
     measured in PennyLane.
@@ -812,7 +827,7 @@ def molecular_hamiltonian(
     args=None,
     load_data=False,
     convert_tol=1e012,
-):  # pylint:disable=too-many-arguments
+):  # pylint:disable=too-many-arguments, too-many-statements
     r"""Generate the qubit Hamiltonian of a molecule.
 
     This function drives the construction of the second-quantized electronic Hamiltonian
@@ -850,8 +865,9 @@ def molecular_hamiltonian(
         basis (str): atomic basis set used to represent the molecular orbitals
         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, or ``method="pyscf"``
-            to use the OpenFermion-PySCF plugin (this requires ``openfermionpyscf`` to be installed).
+            to specify the built-in differentiable Hartree-Fock solver, ``method="pyscf"`` to use
+            the PySCF package (requires ``pyscf`` to be installed), or ``method="openfermion"`` to
+            use the OpenFermion-PySCF plugin (this requires ``openfermionpyscf`` to be installed).
         active_electrons (int): Number of active electrons. If not specified, all electrons
             are considered to be active.
         active_orbitals (int): Number of active orbitals. If not specified, all orbitals
@@ -899,8 +915,8 @@ def molecular_hamiltonian(
     + (0.176276408043196) [Z2 Z3]
     """
 
-    if method not in ["dhf", "pyscf"]:
-        raise ValueError("Only 'dhf' and 'pyscf' backends are supported.")
+    if method not in ["dhf", "pyscf", "openfermion"]:
+        raise ValueError("Only 'dhf', 'pyscf' and 'openfermion' backends are supported.")
 
     if len(coordinates) == len(symbols) * 3:
         geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3))
@@ -909,10 +925,13 @@ def molecular_hamiltonian(
         geometry_dhf = qml.numpy.array(coordinates)
         geometry_hf = coordinates.flatten()
 
+    wires_map = None
+
+    if wires:
+        wires_new = qml.qchem.convert._process_wires(wires)
+        wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels)))
+
     if method == "dhf":
-        if wires:
-            wires_new = qml.qchem.convert._process_wires(wires)
-            wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels)))
 
         if mapping != "jordan_wigner":
             raise ValueError(
@@ -964,9 +983,25 @@ def molecular_hamiltonian(
             h = qml.map_wires(h, wires_map)
         return h, 2 * len(active)
 
+    if method == "pyscf" and mapping.strip().lower() == "jordan_wigner":
+        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)
+
+        if active_new_opmath():
+            h_pl = qml.jordan_wigner(hf, wire_map=wires_map, tol=1.0e-10).simplify()
+
+        else:
+            h_pl = qml.jordan_wigner(hf, ps=True, wire_map=wires_map, tol=1.0e-10).hamiltonian()
+            h_pl = simplify(h_pl)
+
+        return h_pl, len(h_pl.wires)
+
     openfermion, _ = _import_of()
 
-    hf_file = meanfield(symbols, geometry_hf, name, charge, mult, basis, method, outpath)
+    hf_file = meanfield(symbols, geometry_hf, name, charge, mult, basis, "pyscf", outpath)
 
     molecule = openfermion.MolecularData(filename=hf_file)
 
@@ -978,4 +1013,73 @@ def molecular_hamiltonian(
 
     h_pl = qml.qchem.convert.import_operator(h_of, wires=wires, tol=convert_tol)
 
-    return h_pl, qubits
+    return h_pl, len(h_pl.wires)
+
+
+def _pyscf_integrals(
+    symbols,
+    coordinates,
+    charge=0,
+    mult=1,
+    basis="sto-3g",
+    active_electrons=None,
+    active_orbitals=None,
+):
+    r"""Compute pyscf integrals."""
+
+    pyscf = _import_pyscf()
+
+    geometry = [
+        [symbol, tuple(np.array(coordinates)[3 * i : 3 * i + 3] * bohr_angs)]
+        for i, symbol in enumerate(symbols)
+    ]
+
+    # create the Mole object
+    mol = pyscf.gto.Mole()
+    mol.atom = geometry
+    mol.basis = basis
+    mol.charge = charge
+    mol.spin = mult - 1
+    mol.verbose = 0
+
+    # initialize the molecule
+    mol_pyscf = mol.build()
+
+    # run HF calculations
+    if mult == 1:
+        molhf = pyscf.scf.RHF(mol_pyscf)
+    else:
+        molhf = pyscf.scf.ROHF(mol_pyscf)
+    _ = molhf.kernel()
+
+    # compute atomic and molecular orbitals
+    one_ao = mol_pyscf.intor_symmetric("int1e_kin") + mol_pyscf.intor_symmetric("int1e_nuc")
+    two_ao = mol_pyscf.intor("int2e_sph")
+    one_mo = np.einsum("pi,pq,qj->ij", molhf.mo_coeff, one_ao, molhf.mo_coeff, optimize=True)
+    two_mo = pyscf.ao2mo.incore.full(two_ao, molhf.mo_coeff)
+
+    core_constant = np.array([molhf.energy_nuc()])
+
+    # convert the two-electron integral tensor to the physicists’ notation
+    two_mo = np.swapaxes(two_mo, 1, 3)
+
+    # define the active space and recompute the integrals
+    core, active = qml.qchem.active_space(
+        mol.nelectron, mol.nao, mult, active_electrons, active_orbitals
+    )
+
+    if core and active:
+        for i in core:
+            core_constant = core_constant + 2 * one_mo[i][i]
+            for j in core:
+                core_constant = core_constant + 2 * two_mo[i][j][j][i] - two_mo[i][j][i][j]
+
+        for p in active:
+            for q in active:
+                for i in core:
+                    one_mo[p, q] = one_mo[p, q] + (2 * two_mo[i][p][q][i] - two_mo[i][p][i][q])
+
+        one_mo = one_mo[qml.math.ix_(active, active)]
+        two_mo = two_mo[qml.math.ix_(active, active, active, active)]
+
+    return core_constant, one_mo, two_mo


Test patch:
diff --git a/tests/qchem/of_tests/test_molecular_hamiltonian.py b/tests/qchem/of_tests/test_molecular_hamiltonian.py
index 18c42f2a0..460f6aaa8 100644
--- a/tests/qchem/of_tests/test_molecular_hamiltonian.py
+++ b/tests/qchem/of_tests/test_molecular_hamiltonian.py
@@ -14,7 +14,7 @@
 """
 Unit tests for molecular Hamiltonians.
 """
-# pylint: disable=too-many-arguments
+# pylint: disable=too-many-arguments, protected-access
 import pytest
 
 from pennylane import Identity, PauliX, PauliY, PauliZ
@@ -388,7 +388,9 @@ def test_diff_hamiltonian_error(symbols, geometry):
     with pytest.raises(ValueError, match="Only 'jordan_wigner' mapping is supported"):
         qchem.molecular_hamiltonian(symbols, geometry, method="dhf", mapping="bravyi_kitaev")
 
-    with pytest.raises(ValueError, match="Only 'dhf' and 'pyscf' backends are supported"):
+    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"):
@@ -440,3 +442,37 @@ def test_real_hamiltonian(symbols, geometry, method, args, tmpdir, op_arithmetic
 
     else:
         assert np.isrealobj(hamiltonian.coeffs)
+
+
+@pytest.mark.parametrize(
+    ("symbols", "geometry", "core_ref", "one_ref", "two_ref"),
+    [
+        (
+            ["H", "H"],
+            np.array([0.0, 0.0, 0.0, 0.0, 0.0, 2.0]),
+            np.array([0.5]),
+            np.array([[-1.08269537e00, 1.88626892e-13], [1.88848936e-13, -6.04947784e-01]]),
+            np.array(
+                [
+                    [
+                        [[6.16219836e-01, -1.93289829e-13], [-1.93373095e-13, 2.00522469e-01]],
+                        [[-1.93345340e-13, 2.00522469e-01], [6.13198399e-01, -1.86684002e-13]],
+                    ],
+                    [
+                        [[-1.93289829e-13, 6.13198399e-01], [2.00522469e-01, -1.86572979e-13]],
+                        [[2.00522469e-01, -1.86961557e-13], [-1.86684002e-13, 6.43874664e-01]],
+                    ],
+                ]
+            ),
+        ),
+    ],
+)
+@pytest.mark.usefixtures("skip_if_no_openfermion_support")
+def test_pyscf_integrals(symbols, geometry, core_ref, one_ref, two_ref):
+    r"""Test that _pyscf_integrals returns correct integrals."""
+
+    core, one, two = qchem.openfermion_obs._pyscf_integrals(symbols, geometry)
+
+    assert np.allclose(core, core_ref)
+    assert np.allclose(one, one_ref)
+    assert np.allclose(two, two_ref)
diff --git a/tests/qchem/of_tests/test_observable_of.py b/tests/qchem/of_tests/test_observable_of.py
index 30e240391..a74d37d4d 100644
--- a/tests/qchem/of_tests/test_observable_of.py
+++ b/tests/qchem/of_tests/test_observable_of.py
@@ -180,3 +180,14 @@ def test_import_of(monkeypatch):
 
         with pytest.raises(ImportError, match="This feature requires openfermion"):
             qml.qchem.openfermion_obs._import_of()
+
+
+def test_import_pyscf(monkeypatch):
+    """Test if an ImportError is raised by _import_pyscf function."""
+    # pylint: disable=protected-access
+
+    with monkeypatch.context() as m:
+        m.setitem(sys.modules, "pyscf", None)
+
+        with pytest.raises(ImportError, match="This feature requires pyscf"):
+            qml.qchem.openfermion_obs._import_pyscf()