Operator Arithmetic support to molecular_hamiltonian() in qchem (#4159

)

* Added warning and op math support

* docs

* more tests

* finished testing

* Apply suggestions from code review

* Apply suggestions from code review

Co-authored-by: Utkarsh <utkarshazad98@gmail.com>

* address code review comments

---------

Co-authored-by: Utkarsh <utkarshazad98@gmail.com>

doc/development/deprecations.rst

@@ -6,6 +6,23 @@ Deprecations
 Pending deprecations
 --------------------
 
+* The `grouping_type` and `grouping_method` arguments of `qchem.molecular_hamiltonian()` are deprecated.
+
+ - Deprecated in v0.31
+ - Will be removed in v0.32
+
+ Instead, simply construct a new instance of ``Hamiltonian`` with the grouping specified:
+
+ .. code-block:: python
+
+ H, qubits = molecular_hamiltonian(symbols, coordinates)
+ grouped_h = qml.Hamiltonian(
+ H.coeffs,
+ H.ops,
+ grouping_type=grouping_type,
+ method=grouping_method,
+ )
+
 * ``zyz_decomposition`` and ``xyx_decomposition`` are deprecated, use ``one_qubit_decomposition`` with a rotations
  keyword instead.
 

doc/releases/changelog-dev.md

@@ -169,6 +169,9 @@
 * The construction of the pauli representation for the `Sum` class is now faster.
  [(#4142)](https://github.com/PennyLaneAI/pennylane/pull/4142)
 
+* `qchem.molecular_hamiltonian()` will now return an arithmetic operator if `enable_new_opmath()` is active.
+ [(#4159)](https://github.com/PennyLaneAI/pennylane/pull/4159)
+
 * `qchem.qubit_observable()` will now return an arithmetic operator if `enable_new_opmath()` is active. 
  [(#4138)](https://github.com/PennyLaneAI/pennylane/pull/4138)
 

pennylane/qchem/openfermion_obs.py

@@ -16,10 +16,12 @@
 # pylint: disable=too-many-arguments, too-few-public-methods, too-many-branches, unused-variable
 # pylint: disable=consider-using-generator, protected-access
 import os
+import warnings
 
 import numpy as np
 
 import pennylane as qml
+from pennylane.operation import active_new_opmath
 
 # Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0)
 bohr_angs = 0.529177210903
@@ -837,6 +839,13 @@ def molecular_hamiltonian(
 
  |
 
+ .. warning::
+
+ If this function is called with a :code:`grouping_type` while :code:`enable_new_opmath()` is active, then
+ an arithmetic operator will be returned and the grouping arguments will be ignored.
+
+ |
+
  Args:
  symbols (list[str]): symbols of the atomic species in the molecule
  coordinates (array[float]): atomic positions in Cartesian coordinates.
@@ -903,6 +912,22 @@ def molecular_hamiltonian(
  + (0.176276408043196) [Z2 Z3]
  """
 
+ if grouping_type is not None:
+ if active_new_opmath():
+ warnings.warn(
+ "The 'grouping_type' and 'grouping_method' arguments are ignored when "
+ "'active_new_opmath() is True. Please disable it (via 'disable_new_opmath()') "
+ "in order to utilize the grouping functionality.",
+ UserWarning,
+ )
+
+ warnings.warn(
+ "The grouping functionality in molecular_hamiltonian is deprecated. "
+ "Please manually compute the groupings via: "
+ "qml.Hamiltonian(molecular_h.coeffs, molecular_h.ops, grouping_type=grouping_type, method=grouping_method)",
+ UserWarning,
+ )
+
  if method not in ["dhf", "pyscf"]:
  raise ValueError("Only 'dhf' and 'pyscf' backends are supported.")
 
@@ -942,14 +967,28 @@ def molecular_hamiltonian(
  core, active = qml.qchem.active_space(
  mol.n_electrons, mol.n_orbitals, mult, active_electrons, active_orbitals
  )
- if args is None:
- h = qml.qchem.diff_hamiltonian(mol, core=core, active=active)()
- coeffs = qml.numpy.real(h.coeffs, requires_grad=False)
- h = qml.Hamiltonian(coeffs, h.ops, grouping_type=grouping_type, method=grouping_method)
+
+ requires_grad = args is not None
+ h = (
+ qml.qchem.diff_hamiltonian(mol, core=core, active=active)(*args)
+ if requires_grad
+ else qml.qchem.diff_hamiltonian(mol, core=core, active=active)()
+ )
+
+ if active_new_opmath():
+ h_as_ps = qml.pauli.pauli_sentence(h)
+ coeffs = qml.numpy.real(list(h_as_ps.values()), requires_grad=requires_grad)
+
+ h_as_ps = qml.pauli.PauliSentence(dict(zip(h_as_ps.keys(), coeffs)))
+ h = (
+ qml.s_prod(0, qml.Identity(h.wires[0]))
+ if len(h_as_ps) == 0
+ else h_as_ps.operation()
+ )
  else:
- h = qml.qchem.diff_hamiltonian(mol, core=core, active=active)(*args)
- coeffs = qml.numpy.real(h.coeffs)
+ coeffs = qml.numpy.real(h.coeffs, requires_grad=requires_grad)
  h = qml.Hamiltonian(coeffs, h.ops, grouping_type=grouping_type, method=grouping_method)
+
  if wires:
  h = qml.map_wires(h, wires_map)
  return h, 2 * len(active)
@@ -968,7 +1007,4 @@ def molecular_hamiltonian(
 
  h_pl = qml.qchem.convert.import_operator(h_of, wires=wires, tol=convert_tol)
 
- return (
- qml.Hamiltonian(h_pl.coeffs, h_pl.ops, grouping_type=grouping_type, method=grouping_method),
- qubits,
- )
+ return h_pl, qubits

tests/qchem/of_tests/test_molecular_hamiltonian.py

@@ -1,11 +1,15 @@
 import sys
+import warnings
 
 import pytest
 
 from pennylane import Identity, PauliX, PauliY, PauliZ
 from pennylane import numpy as np
 from pennylane import qchem
 from pennylane.ops import Hamiltonian
+from pennylane.ops.functions import dot
+from pennylane.pauli import pauli_sentence
+from pennylane.operation import enable_new_opmath, disable_new_opmath
 
 # TODO: Bring pytest skip to relevant tests.
 openfermion = pytest.importorskip("openfermion")
@@ -42,6 +46,7 @@
 )
 
 
+@pytest.mark.parametrize("op_arithmetic", [False, True])
 @pytest.mark.parametrize(
  (
  "charge",
@@ -68,28 +73,44 @@ def test_building_hamiltonian(
  mapping,
  grouping_type,
  tmpdir,
+ op_arithmetic,
 ):
  r"""Test that the generated Hamiltonian `built_hamiltonian` 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"""
- built_hamiltonian, qubits = qchem.molecular_hamiltonian(
- symbols,
- coordinates,
- charge=charge,
- mult=mult,
- method=package,
- active_electrons=nact_els,
- active_orbitals=nact_orbs,
- mapping=mapping,
- grouping_type=grouping_type,
- outpath=tmpdir.strpath,
- )
 
- assert isinstance(built_hamiltonian, Hamiltonian)
+ args = (symbols, coordinates)
+ kwargs = {
+ "charge": charge,
+ "mult": mult,
+ "method": package,
+ "active_electrons": nact_els,
+ "active_orbitals": nact_orbs,
+ "mapping": mapping,
+ "grouping_type": grouping_type,
+ "outpath": tmpdir.strpath,
+ }
+ if op_arithmetic:
+ enable_new_opmath()
+
+ if grouping_type is not None:
+ with pytest.warns(
+ UserWarning, match="The grouping functionality in molecular_hamiltonian is deprecated."
+ ):
+ built_hamiltonian, qubits = qchem.molecular_hamiltonian(*args, **kwargs)
+ else:
+ built_hamiltonian, qubits = qchem.molecular_hamiltonian(*args, **kwargs)
+
+ if op_arithmetic:
+ disable_new_opmath()
+ assert not isinstance(built_hamiltonian, Hamiltonian)
+ else:
+ assert isinstance(built_hamiltonian, Hamiltonian)
  assert qubits == 2 * nact_orbs
 
 
+@pytest.mark.parametrize("op_arithmetic", [False, True])
 @pytest.mark.parametrize(
  ("symbols", "geometry", "h_ref_data"),
  [
@@ -189,32 +210,56 @@ def test_building_hamiltonian(
  ),
  ],
 )
-def test_differentiable_hamiltonian(symbols, geometry, h_ref_data):
+def test_differentiable_hamiltonian(symbols, geometry, h_ref_data, op_arithmetic):
  r"""Test that molecular_hamiltonian returns the correct Hamiltonian with the differentiable
  backend."""
+ if op_arithmetic:
+ enable_new_opmath()
+
  geometry.requires_grad = True
  args = [geometry.reshape(2, 3)]
  h_args = qchem.molecular_hamiltonian(symbols, geometry, method="dhf", args=args)[0]
 
  geometry.requires_grad = False
  h_noargs = qchem.molecular_hamiltonian(symbols, geometry, method="dhf", grouping_type="qwc")[0]
 
- h_ref = Hamiltonian(h_ref_data[0], h_ref_data[1])
-
- assert np.allclose(np.sort(h_args.coeffs), np.sort(h_ref.coeffs))
- assert Hamiltonian(np.ones(len(h_args.coeffs)), h_args.ops).compare(
- Hamiltonian(np.ones(len(h_ref.coeffs)), h_ref.ops)
+ h_ref = (
+ dot(h_ref_data[0], h_ref_data[1], pauli=True)
+ if op_arithmetic
+ else Hamiltonian(h_ref_data[0], h_ref_data[1])
  )
 
- assert np.allclose(np.sort(h_noargs.coeffs), np.sort(h_ref.coeffs))
- assert Hamiltonian(np.ones(len(h_noargs.coeffs)), h_noargs.ops).compare(
- Hamiltonian(np.ones(len(h_ref.coeffs)), h_ref.ops)
- )
+ if op_arithmetic:
+ disable_new_opmath()
+ h_args_ps = pauli_sentence(h_args)
+ h_noargs_ps = pauli_sentence(h_noargs)
+ h_ref_sorted_coeffs = np.sort(list(h_ref.values()))
+
+ assert set(h_args_ps) == set(h_ref)
+ assert set(h_noargs_ps) == set(h_ref)
+
+ assert np.allclose(np.sort(list(h_args_ps.values())), h_ref_sorted_coeffs)
+ assert np.allclose(np.sort(list(h_noargs_ps.values())), h_ref_sorted_coeffs)
+
+ assert all([val.requires_grad is True for val in h_args_ps.values()])
+ assert all([val.requires_grad is False for val in h_noargs_ps.values()])
 
- assert h_args.coeffs.requires_grad == True
- assert h_noargs.coeffs.requires_grad == False
+ else:
+ assert np.allclose(np.sort(h_args.coeffs), np.sort(h_ref.coeffs))
+ assert Hamiltonian(np.ones(len(h_args.coeffs)), h_args.ops).compare(
+ Hamiltonian(np.ones(len(h_ref.coeffs)), h_ref.ops)
+ )
 
+ assert np.allclose(np.sort(h_noargs.coeffs), np.sort(h_ref.coeffs))
+ assert Hamiltonian(np.ones(len(h_noargs.coeffs)), h_noargs.ops).compare(
+ Hamiltonian(np.ones(len(h_ref.coeffs)), h_ref.ops)
+ )
 
+ assert h_args.coeffs.requires_grad == True
+ assert h_noargs.coeffs.requires_grad == False
+
+
+@pytest.mark.parametrize("op_arithmetic", [False, True])
 @pytest.mark.parametrize(
  ("symbols", "geometry", "method", "wiremap", "grouping"),
  [
@@ -234,8 +279,13 @@ def test_differentiable_hamiltonian(symbols, geometry, h_ref_data):
  ),
  ],
 )
-def test_custom_wiremap_hamiltonian_pyscf(symbols, geometry, method, wiremap, grouping, tmpdir):
+def test_custom_wiremap_hamiltonian_pyscf(
+ symbols, geometry, method, wiremap, grouping, tmpdir, op_arithmetic
+):
  r"""Test that the generated Hamiltonian has the correct wire labels given by a custom wiremap."""
+ if op_arithmetic:
+ enable_new_opmath()
+
  hamiltonian, qubits = qchem.molecular_hamiltonian(
  symbols=symbols,
  coordinates=geometry,
@@ -247,7 +297,11 @@ def test_custom_wiremap_hamiltonian_pyscf(symbols, geometry, method, wiremap, gr
 
  assert set(hamiltonian.wires) == set(wiremap)
 
+ if op_arithmetic:
+ disable_new_opmath()
 
+
+@pytest.mark.parametrize("op_arithmetic", [False, True])
 @pytest.mark.parametrize(
  ("symbols", "geometry", "wiremap", "args", "grouping"),
  [
@@ -267,8 +321,13 @@ def test_custom_wiremap_hamiltonian_pyscf(symbols, geometry, method, wiremap, gr
  ),
  ],
 )
-def test_custom_wiremap_hamiltonian_dhf(symbols, geometry, wiremap, args, grouping, tmpdir):
+def test_custom_wiremap_hamiltonian_dhf(
+ symbols, geometry, wiremap, args, grouping, tmpdir, op_arithmetic
+):
  r"""Test that the generated Hamiltonian has the correct wire labels given by a custom wiremap."""
+ if op_arithmetic:
+ enable_new_opmath()
+
  hamiltonian, qubits = qchem.molecular_hamiltonian(
  symbols=symbols,
  coordinates=geometry,
@@ -280,6 +339,9 @@ def test_custom_wiremap_hamiltonian_dhf(symbols, geometry, wiremap, args, groupi
 
  assert list(hamiltonian.wires) == list(wiremap)
 
+ if op_arithmetic:
+ disable_new_opmath()
+
 
 file_content = """\
 2
@@ -326,6 +388,32 @@ def test_diff_hamiltonian_error(symbols, geometry):
  qchem.molecular_hamiltonian(symbols, geometry, mult=3)[0]
 
 
+@pytest.mark.parametrize(
+ ("symbols", "geometry"),
+ [
+ (
+ ["H", "H"],
+ np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]),
+ ),
+ ],
+)
+def test_hamiltonian_warnings(symbols, geometry):
+ r"""Test that molecular_hamiltonian raises a warning when using grouping."""
+
+ with pytest.warns(
+ UserWarning, match="The grouping functionality in molecular_hamiltonian is deprecated. "
+ ):
+ qchem.molecular_hamiltonian(symbols, geometry, grouping_type="qwc")
+
+ with pytest.warns(
+ UserWarning, match="The 'grouping_type' and 'grouping_method' arguments are ignored when "
+ ):
+ enable_new_opmath()
+ qchem.molecular_hamiltonian(symbols, geometry, grouping_type="qwc")
+ disable_new_opmath()
+
+
+@pytest.mark.parametrize("op_arithmetic", [False, True])
 @pytest.mark.parametrize(
  ("symbols", "geometry", "method", "args"),
  [
@@ -349,8 +437,11 @@ def test_diff_hamiltonian_error(symbols, geometry):
  ),
  ],
 )
-def test_real_hamiltonian(symbols, geometry, method, args, tmpdir):
+def test_real_hamiltonian(symbols, geometry, method, args, tmpdir, op_arithmetic):
  r"""Test that the generated Hamiltonian has real coefficients."""
+ if op_arithmetic:
+ enable_new_opmath()
+
  hamiltonian, qubits = qchem.molecular_hamiltonian(
  symbols=symbols,
  coordinates=geometry,
@@ -359,4 +450,10 @@ def test_real_hamiltonian(symbols, geometry, method, args, tmpdir):
  outpath=tmpdir.strpath,
  )
 
- assert np.isrealobj(hamiltonian.coeffs)
+ if op_arithmetic:
+ disable_new_opmath()
+ h_as_ps = pauli_sentence(hamiltonian)
+ assert np.isrealobj(np.array(h_as_ps.values()))
+
+ else:
+ assert np.isrealobj(hamiltonian.coeffs)