From 283c7afffdc707b7a27a4a540516c08320f40810 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Fri, 3 May 2024 09:32:27 -0400 Subject: [PATCH 01/28] New branch for qchem features --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 1f049ef7616..3c3f533c80a 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -890,7 +890,7 @@ def molecular_hamiltonian( 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 **Example** From 456cd573dc8a647738fde303a2e793a268411ee3 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Fri, 17 May 2024 10:50:37 -0400 Subject: [PATCH 02/28] Make Molecule the central object for qchem functions (#5571) **Context:** **Description of the Change:** Molecule is the central object for molecular_hamiltonian function **Benefits:** **Possible Drawbacks:** **Related GitHub Issues:** --------- Co-authored-by: soranjh Co-authored-by: Utkarsh Co-authored-by: Austin Huang <65315367+austingmhuang@users.noreply.github.com> --- doc/code/qml_qchem.rst | 15 +- pennylane/qchem/molecule.py | 3 +- pennylane/qchem/openfermion_obs.py | 208 ++++++-- requirements-ci.txt | 2 +- .../of_tests/test_molecular_hamiltonian.py | 477 +++++++++++++++--- 5 files changed, 594 insertions(+), 111 deletions(-) diff --git a/doc/code/qml_qchem.rst b/doc/code/qml_qchem.rst index 68103da2831..c4d1b25e9a0 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 @@ -183,14 +184,8 @@ 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' - ) + molecule = qml.qchem.Molecule(symbols, geometry, charge=0, mult=1, basis='sto-3g') + hamiltonian, qubits = qml.qchem.molecular_hamiltonian(molecule, method='pyscf') The non-differentiable backend requires the ``OpenFermion-PySCF`` plugin to be installed by the user with diff --git a/pennylane/qchem/molecule.py b/pennylane/qchem/molecule.py index 94a2fb101b5..24e67d43da0 100644 --- a/pennylane/qchem/molecule.py +++ b/pennylane/qchem/molecule.py @@ -69,6 +69,7 @@ def __init__( charge=0, mult=1, basis_name="sto-3g", + name="molecule", load_data=False, l=None, alpha=None, @@ -99,7 +100,7 @@ def __init__( self.charge = charge self.mult = mult self.basis_name = basis_name.lower() - + self.name = name self.n_basis, self.basis_data = mol_basis_data(self.basis_name, self.symbols, load_data) self.nuclear_charges = [atomic_numbers[s] for s in self.symbols] diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 3c3f533c80a..69716d7903d 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -16,12 +16,14 @@ # 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 # Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0) bohr_angs = 0.529177210903 @@ -811,26 +813,10 @@ def decompose(hf_file, mapping="jordan_wigner", core=None, active=None): 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, load_data=False, 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 +840,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 +857,28 @@ 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 `_ 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 and the number of qubits + .. warning:: + Use of ``qml.qchem.molecular_hamiltonian`` with symbols and geometry arguments is being deprecated. + Instead, please use the method with ``qml.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,8 +897,164 @@ 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, + load_data=False, + 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, + load_data, + convert_tol, + ) + + +@_molecular_hamiltonian_dispatch.register(list) +def _( + 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, +): + 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.") diff --git a/requirements-ci.txt b/requirements-ci.txt index dafa106ffe4..9f7c3b39c77 100644 --- a/requirements-ci.txt +++ b/requirements-ci.txt @@ -8,7 +8,7 @@ autograd toml appdirs semantic_version -autoray +autoray>=0.6.1,<0.6.10 matplotlib requests tomli # Drop once minimum Python version is 3.11 diff --git a/tests/qchem/of_tests/test_molecular_hamiltonian.py b/tests/qchem/of_tests/test_molecular_hamiltonian.py index fc449b3a2c7..8d8d6850875 100644 --- a/tests/qchem/of_tests/test_molecular_hamiltonian.py +++ b/tests/qchem/of_tests/test_molecular_hamiltonian.py @@ -18,7 +18,7 @@ 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 @@ -105,6 +105,53 @@ def test_building_hamiltonian( assert qubits == 2 * nact_orbs +@pytest.mark.parametrize( + ( + "charge", + "mult", + "package", + "nact_els", + "nact_orbs", + "mapping", + ), + [ + (0, 1, "pyscf", 2, 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", "h_ref_data"), [ @@ -137,21 +184,21 @@ 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), ], ), ), @@ -184,21 +231,21 @@ 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), ], ), ), @@ -240,28 +287,159 @@ 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", "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]), + # 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( + [ + 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]]), + # 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( + [ + 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, 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)[0] + + geometry.requires_grad = False + molecule = qchem.Molecule(symbols, geometry) + h_noargs = qchem.molecular_hamiltonian(molecule, method="dhf")[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 +451,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 +517,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,18 +580,29 @@ 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.""" + symbols = ["H", "H"] + geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]) + with pytest.raises(ValueError, match="Only 'jordan_wigner' mapping is supported"): qchem.molecular_hamiltonian(symbols, geometry, method="dhf", mapping="bravyi_kitaev") @@ -363,33 +615,46 @@ def test_diff_hamiltonian_error(symbols, geometry): qchem.molecular_hamiltonian(symbols, geometry, mult=3) +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]) + + with pytest.raises(ValueError, match="Only 'jordan_wigner' mapping is supported"): + qchem.molecular_hamiltonian(symbols, geometry, method="dhf", mapping="bravyi_kitaev") + + molecule = qchem.Molecule(symbols, geometry) + with pytest.raises( + ValueError, match="Only 'dhf', 'pyscf' and 'openfermion' backends are supported" + ): + qchem.molecular_hamiltonian(molecule, method="psi4") + + @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 +666,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 +733,52 @@ 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") From 3e3770923be61e212e23f116ae51d997d4c67f09 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Fri, 24 May 2024 10:06:58 -0400 Subject: [PATCH 03/28] Support Angstrom units with Molecule class and molecular_hamiltonian function (#5694) **Context:** Support Angstrom units along with Bohr **Description of the Change:** Molecule class and molecular_hamiltonian function now works with both Bohr and Angstrom units. **Benefits:** **Possible Drawbacks:** **Related GitHub Issues:** --------- Co-authored-by: soranjh Co-authored-by: Utkarsh 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/releases/changelog-dev.md | 4 + pennylane/qchem/molecule.py | 15 ++++ pennylane/qchem/openfermion_obs.py | 11 +++ requirements-ci.txt | 4 +- .../of_tests/test_molecular_hamiltonian.py | 74 +++++++++++++++++++ tests/qchem/test_molecule.py | 9 +++ 6 files changed, 115 insertions(+), 2 deletions(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 9d88744ff53..a0418b58e1b 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -80,6 +80,10 @@ `qml.devices.Device`, which follows the new device API. [(#5581)](https://github.com/PennyLaneAI/pennylane/pull/5581) +* `qml.qchem.Molecule` and qchem functions that take `Molecule` object as an argument now work with coordinates + provided in units Angstrom along with Bohr. + [(#5694)](https://github.com/PennyLaneAI/pennylane/pull/5694) + * The `dtype` for `eigvals` of `X`, `Y`, `Z` and `Hadamard` is changed from `int` to `float`, making them consistent with the other observables. The `dtype` of the returned values when sampling these observables (e.g. `qml.sample(X(0))`) is also changed to `float`. diff --git a/pennylane/qchem/molecule.py b/pennylane/qchem/molecule.py index 24e67d43da0..7f6877a0fbb 100644 --- a/pennylane/qchem/molecule.py +++ b/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. @@ -51,6 +54,7 @@ class Molecule: 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** @@ -75,6 +79,7 @@ def __init__( alpha=None, coeff=None, normalize=True, + unit="bohr", ): if ( basis_name.lower() @@ -97,12 +102,22 @@ 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.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 diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 69716d7903d..a9252e05505 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -997,6 +997,7 @@ def _( def _( symbols, coordinates, + unit="bohr", name="molecule", charge=0, mult=1, @@ -1013,6 +1014,16 @@ def _( 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, diff --git a/requirements-ci.txt b/requirements-ci.txt index 6ff2c5afc86..932820a266d 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 diff --git a/tests/qchem/of_tests/test_molecular_hamiltonian.py b/tests/qchem/of_tests/test_molecular_hamiltonian.py index 8d8d6850875..355a6b19b6a 100644 --- a/tests/qchem/of_tests/test_molecular_hamiltonian.py +++ b/tests/qchem/of_tests/test_molecular_hamiltonian.py @@ -782,3 +782,77 @@ def test_error_raised_for_missing_molecule_information(): 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( + ("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_molecule.py b/tests/qchem/test_molecule.py index 5191f3445f7..bd01511cbe9 100644 --- a/tests/qchem/test_molecule.py +++ b/tests/qchem/test_molecule.py @@ -353,3 +353,12 @@ def test_load_basis(self, symbols, geometry, basis_name, params_ref): 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") From 8b982a065b129e4b511229e7f70ad69ac89c49fe Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Fri, 24 May 2024 13:32:16 -0400 Subject: [PATCH 04/28] correct openshell errors in qchem (#5655) **Context:** Changes openshell warnings. **Description of the Change:** Changed openshell warnings to work with new molecular_hamiltonian function **Benefits:** molecular_hamiltonian function works with openshell molecules when method is pyscf **Possible Drawbacks:** **Related GitHub Issues:** --------- Co-authored-by: soranjh Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com> Co-authored-by: Austin Huang <65315367+austingmhuang@users.noreply.github.com> Co-authored-by: soranjh <40344468+soranjh@users.noreply.github.com> --- pennylane/qchem/hartree_fock.py | 5 ++++ pennylane/qchem/molecule.py | 9 +----- pennylane/qchem/openfermion_obs.py | 30 ++++++++----------- .../of_tests/test_molecular_hamiltonian.py | 22 ++++++++++++-- tests/qchem/test_hartree_fock.py | 10 +++++++ tests/qchem/test_molecule.py | 11 ------- 6 files changed, 47 insertions(+), 40 deletions(-) diff --git a/pennylane/qchem/hartree_fock.py b/pennylane/qchem/hartree_fock.py index 381183372ae..69750cc8d07 100644 --- a/pennylane/qchem/hartree_fock.py +++ b/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 diff --git a/pennylane/qchem/molecule.py b/pennylane/qchem/molecule.py index 7f6877a0fbb..1dd831bf9d3 100644 --- a/pennylane/qchem/molecule.py +++ b/pennylane/qchem/molecule.py @@ -44,8 +44,7 @@ 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'. load_data (bool): flag to load data from the basis-set-exchange library @@ -122,12 +121,6 @@ def __init__( 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/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index a9252e05505..0a8cfd469f3 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -25,6 +25,8 @@ 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 @@ -594,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;" @@ -1082,17 +1071,22 @@ 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 in ("dhf", "pyscf"): + n_electrons = sum([atomic_numbers[s] for s in symbols]) - charge + + if n_electrons % 2 == 1 or mult != 1: + raise ValueError( + "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 mapping != "jordan_wigner": raise ValueError( "Only 'jordan_wigner' mapping is supported for the differentiable workflow." ) - if 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." - ) + if args is None and isinstance(geometry_dhf, qml.numpy.tensor): geometry_dhf.requires_grad = False mol = qml.qchem.Molecule( diff --git a/tests/qchem/of_tests/test_molecular_hamiltonian.py b/tests/qchem/of_tests/test_molecular_hamiltonian.py index 355a6b19b6a..6da338e43fb 100644 --- a/tests/qchem/of_tests/test_molecular_hamiltonian.py +++ b/tests/qchem/of_tests/test_molecular_hamiltonian.py @@ -65,8 +65,8 @@ ), [ (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"), ], ) @@ -116,6 +116,8 @@ def test_building_hamiltonian( ), [ (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"), ], ) @@ -611,10 +613,24 @@ def test_diff_hamiltonian_error(): ): 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.""" diff --git a/tests/qchem/test_hartree_fock.py b/tests/qchem/test_hartree_fock.py index a0bcf2b5685..2709a528a43 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 bd01511cbe9..b4e797e5687 100644 --- a/tests/qchem/test_molecule.py +++ b/tests/qchem/test_molecule.py @@ -46,17 +46,6 @@ def test_basis_error(self, symbols, geometry): 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"), [ From 793d1025b5cc01242a66ff55992bb8011c36d264 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 27 May 2024 14:34:16 -0400 Subject: [PATCH 05/28] Update mapping support in qchem (#5657) **Context:** added different mappings to molecular_hamiltonian function **Description of the Change:** molecular Hamiltonian can be mapped using different mapping schemes. **Benefits:** **Possible Drawbacks:** **Related GitHub Issues:** --------- Co-authored-by: soranjh --- doc/releases/changelog-dev.md | 4 + pennylane/qchem/hamiltonian.py | 6 +- pennylane/qchem/observable_hf.py | 15 +- pennylane/qchem/openfermion_obs.py | 48 +- tests/qchem/of_tests/test_decompose.py | 14 - .../of_tests/test_molecular_hamiltonian.py | 498 +++++++++++++++++- 6 files changed, 530 insertions(+), 55 deletions(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index a0418b58e1b..0ac4a758940 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -103,6 +103,9 @@ * Empty initialization of `PauliVSpace` is permitted. [(#5675)](https://github.com/PennyLaneAI/pennylane/pull/5675) +* `qml.qchem.molecular_hamiltonian function now works with parity and Bravyi-Kitaev mappings. + [(#5657)](https://github.com/PennyLaneAI/pennylane/pull/5657/) +

Community contributions 🥳

* Implemented kwargs (`check_interface`, `check_trainability`, `rtol` and `atol`) support in `qml.equal` for the operators `Pow`, `Adjoint`, `Exp`, and `SProd`. @@ -189,6 +192,7 @@ Lillian M. A. Frederiksen, Ahmed Darwish, Gabriel Bottrill, Isaac De Vlugt, +Diksha Dhawan, Pietropaolo Frisoni, Soran Jahangiri, Korbinian Kottmann, diff --git a/pennylane/qchem/hamiltonian.py b/pennylane/qchem/hamiltonian.py index 629bc229ae5..f0d6696a475 100644 --- a/pennylane/qchem/hamiltonian.py +++ b/pennylane/qchem/hamiltonian.py @@ -181,7 +181,7 @@ 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: @@ -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 _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 diff --git a/pennylane/qchem/observable_hf.py b/pennylane/qchem/observable_hf.py index ffc5326d256..bc18dc18288 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 0a8cfd469f3..ddb26670043 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -789,13 +789,10 @@ 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) @@ -1058,6 +1055,12 @@ def _molecular_hamiltonian( 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 @@ -1082,11 +1085,6 @@ def _molecular_hamiltonian( if method == "dhf": - if mapping != "jordan_wigner": - raise ValueError( - "Only 'jordan_wigner' mapping is supported for the differentiable workflow." - ) - if args is None and isinstance(geometry_dhf, qml.numpy.tensor): geometry_dhf.requires_grad = False mol = qml.qchem.Molecule( @@ -1105,9 +1103,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(): @@ -1128,19 +1126,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/tests/qchem/of_tests/test_decompose.py b/tests/qchem/of_tests/test_decompose.py index ccb4f431439..88adde1f01c 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 6da338e43fb..f412f081436 100644 --- a/tests/qchem/of_tests/test_molecular_hamiltonian.py +++ b/tests/qchem/of_tests/test_molecular_hamiltonian.py @@ -155,16 +155,16 @@ def test_building_hamiltonian_molecule_class( @pytest.mark.parametrize( - ("symbols", "geometry", "h_ref_data"), + ("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( [ @@ -204,14 +204,109 @@ def test_building_hamiltonian_molecule_class( ], ), ), + ( + ["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( [ @@ -254,16 +349,18 @@ def test_building_hamiltonian_molecule_class( ], ) @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 @@ -290,16 +387,16 @@ def test_differentiable_hamiltonian(symbols, geometry, h_ref_data): @pytest.mark.parametrize( - ("symbols", "geometry", "h_ref_data"), + ("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( [ @@ -339,14 +436,109 @@ def test_differentiable_hamiltonian(symbols, geometry, h_ref_data): ], ), ), + ( + ["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( [ @@ -389,18 +581,18 @@ def test_differentiable_hamiltonian(symbols, geometry, h_ref_data): ], ) @pytest.mark.usefixtures("use_legacy_and_new_opmath") -def test_differentiable_hamiltonian_molecule_class(symbols, geometry, h_ref_data): +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)[0] + 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")[0] + h_noargs = qchem.molecular_hamiltonian(molecule, method="dhf", mapping=mapping)[0] ops = [ qml.operation.Tensor(*op) if isinstance(op, qml.ops.Prod) else op @@ -605,9 +797,6 @@ def test_diff_hamiltonian_error(): symbols = ["H", "H"] geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]) - 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', 'pyscf' and 'openfermion' backends are supported" ): @@ -637,15 +826,15 @@ def test_diff_hamiltonian_error_molecule_class(): symbols = ["H", "H"] geometry = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.0]) - with pytest.raises(ValueError, match="Only 'jordan_wigner' mapping is supported"): - qchem.molecular_hamiltonian(symbols, geometry, method="dhf", mapping="bravyi_kitaev") - 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( ("method", "args"), @@ -800,6 +989,279 @@ def test_error_raised_for_missing_molecule_information(): 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"), [ From b49d7df33085da09f66f515c4c5d124cef3ae813 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Wed, 29 May 2024 11:51:46 -0400 Subject: [PATCH 06/28] Added molecular_dipole function --- pennylane/qchem/__init__.py | 1 + pennylane/qchem/dipole.py | 6 +- pennylane/qchem/molecule.py | 1 + pennylane/qchem/openfermion_obs.py | 140 +++++++- tests/qchem/of_tests/test_molecular_dipole.py | 326 ++++++++++++++++++ tests/qchem/of_tests/test_observable_of.py | 2 +- 6 files changed, 470 insertions(+), 6 deletions(-) create mode 100644 tests/qchem/of_tests/test_molecular_dipole.py diff --git a/pennylane/qchem/__init__.py b/pennylane/qchem/__init__.py index df724e57602..afe99678a84 100644 --- a/pennylane/qchem/__init__.py +++ b/pennylane/qchem/__init__.py @@ -53,6 +53,7 @@ from .openfermion_obs import ( decompose, dipole_of, + molecular_dipole, meanfield, molecular_hamiltonian, observable, diff --git a/pennylane/qchem/dipole.py b/pennylane/qchem/dipole.py index de17580cf69..d4c36a7505c 100644 --- a/pennylane/qchem/dipole.py +++ b/pennylane/qchem/dipole.py @@ -228,7 +228,7 @@ def _fermionic_dipole(*args): return _fermionic_dipole -def dipole_moment(mol, cutoff=1.0e-16, core=None, active=None): +def dipole_moment(mol, cutoff=1.0e-16, core=None, active=None, mapping="jordan_wigner"): r"""Return a function that computes the qubit dipole moment observable. The dipole operator in the second-quantized form is @@ -277,6 +277,8 @@ def dipole_moment(mol, cutoff=1.0e-16, core=None, active=None): cutoff (float): cutoff value for discarding the negligible dipole moment integrals core (list[int]): indices of the core orbitals active (list[int]): indices of the active orbitals + mapping (str): Specifies the transformation to map the fermionic Hamiltonian to the + Pauli basis. Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. Returns: function: function that computes the qubit dipole moment observable @@ -313,7 +315,7 @@ def _dipole(*args): d = [] d_ferm = fermionic_dipole(mol, cutoff, core, active)(*args) for i in d_ferm: - d.append(qubit_observable(i, cutoff=cutoff)) + d.append(qubit_observable(i, cutoff=cutoff, mapping=mapping)) return d diff --git a/pennylane/qchem/molecule.py b/pennylane/qchem/molecule.py index 1dd831bf9d3..c6b473e6df1 100644 --- a/pennylane/qchem/molecule.py +++ b/pennylane/qchem/molecule.py @@ -106,6 +106,7 @@ def __init__( 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"): diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index ddb26670043..0e41c6feda3 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -113,7 +113,7 @@ def observable(fermion_ops, init_term=0, mapping="jordan_wigner", wires=None): example, this can be used to pass the nuclear-nuclear repulsion energy :math:`V_{nn}` which is typically included in the electronic Hamiltonian of molecules. mapping (str): Specifies the fermion-to-qubit mapping. Input values can - be ``'jordan_wigner'`` or ``'bravyi_kitaev'``. + be ``'jordan_wigner'``, ``'parity'``, or ``'bravyi_kitaev'``. wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator to an observable measurable in a PennyLane ansatz. For types Wires/list/tuple, each item in the iterable represents a wire label @@ -140,10 +140,10 @@ def observable(fermion_ops, init_term=0, mapping="jordan_wigner", wires=None): """ openfermion, _ = _import_of() - if mapping.strip().lower() not in ("jordan_wigner", "bravyi_kitaev"): + if mapping.strip().lower() not in ("jordan_wigner", "parity", "bravyi_kitaev"): raise TypeError( f"The '{mapping}' transformation is not available. \n " - f"Please set 'mapping' to 'jordan_wigner' or 'bravyi_kitaev'." + f"Please set 'mapping' to 'jordan_wigner', 'parity', or 'bravyi_kitaev'." ) # Initialize the FermionOperator @@ -161,6 +161,15 @@ def observable(fermion_ops, init_term=0, mapping="jordan_wigner", wires=None): openfermion.transforms.bravyi_kitaev(mb_obs), wires=wires ) + if mapping == "parity": + qubits = openfermion.count_qubits(mb_obs) + if qubits == 0: + return 0.0 * qml.I(0) + binary_code = openfermion.parity_code(qubits) + return qml.qchem.convert.import_operator( + openfermion.transforms.binary_code_transform(mb_obs, binary_code), wires=wires + ) + return qml.qchem.convert.import_operator( openfermion.transforms.jordan_wigner(mb_obs), wires=wires ) @@ -649,6 +658,131 @@ def dipole_of( return dip +def molecular_dipole( + molecule, + method="dhf", + active_electrons=None, + active_orbitals=None, + mapping="jordan_wigner", + outpath=".", + wires=None, + args=None, + cutoff=1.0e-16, +): # pylint:disable=too-many-arguments, too-many-statements + r"""Generate the dipole moment operator for a molecule in the Pauli basis. + Args: + 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, `or 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 + are considered to be active. + mapping (str): transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian + outpath (str): path to the directory containing output files + wires (Wires, list, tuple, dict): Custom wire mapping for connecting to Pennylane ansatz. + For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label + 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. + args (array[array[float]]): initial values of the differentiable parameters + cutoff (float): Cutoff value for including the matrix elements + :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements + with absolute value less than ``cutoff`` are neglected. + + Returns: + list[pennylane.Hamiltonian]: the qubit observables corresponding to the components + :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator in + atomic units. + + """ + + if method not in ["dhf", "openfermion"]: + raise ValueError("Only 'dhf', 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'." + ) + + symbols = molecule.symbols + coordinates = molecule.coordinates + + if len(coordinates) == len(symbols) * 3: + print("1D coordinates") + geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3)) + geometry_hf = coordinates + print(geometry_dhf, geometry_hf) + elif len(coordinates) == len(symbols): + print("2D coordinates") + geometry_dhf = qml.numpy.array(coordinates) + geometry_hf = coordinates.flatten() + print(geometry_dhf, geometry_hf) + + if molecule.mult != 1: + raise ValueError( + "Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule." + ) + + 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))) + + core, active = qml.qchem.active_space( + molecule.n_electrons, molecule.n_orbitals, molecule.mult, active_electrons, active_orbitals + ) + + if method == "dhf": + + if args is None and isinstance(geometry_dhf, qml.numpy.tensor): + geometry_dhf.requires_grad = False + mol = qml.qchem.Molecule( + symbols, + geometry_dhf, + charge=molecule.charge, + mult=molecule.mult, + basis_name=molecule.basis_name, + load_data=molecule.load_data, + alpha=molecule.alpha, + coeff=molecule.coeff, + ) + + requires_grad = args is not None + dip = ( + qml.qchem.dipole_moment(mol, cutoff=cutoff, core=core, active=active, mapping=mapping)( + *args + ) + if requires_grad + else qml.qchem.dipole_moment( + mol, cutoff=cutoff, core=core, active=active, mapping=mapping + )() + ) + return dip + + dip = qml.qchem.dipole_of( + symbols, + geometry_hf, + molecule.name, + molecule.charge, + molecule.mult, + molecule.basis_name, + package="pyscf", + core=core, + active=active, + mapping=mapping, + cutoff=cutoff, + outpath=outpath, + wires=wires, + ) + + return dip + + def meanfield( symbols, coordinates, diff --git a/tests/qchem/of_tests/test_molecular_dipole.py b/tests/qchem/of_tests/test_molecular_dipole.py new file mode 100644 index 00000000000..99cd4107b79 --- /dev/null +++ b/tests/qchem/of_tests/test_molecular_dipole.py @@ -0,0 +1,326 @@ +# Copyright 2018-2023 Xanadu Quantum Technologies Inc. + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Unit tests for molecular dipole. +""" +# pylint: disable=too-many-arguments, protected-access +import pytest + +import pennylane as qml +from pennylane import I, X, Y, Z +from pennylane import numpy as np +from pennylane import qchem +from pennylane.operation import active_new_opmath + +h2 = ["H", "H"] +x_h2 = np.array([0.0, 0.0, -0.661, 0.0, 0.0, 0.661]) +coeffs_h2 = [] +coeffs_h2.append([0.0]) +coeffs_h2.append([0.0]) +coeffs_h2.append([0.45445016, 0.45445016, 0.45445016, 0.45445016]) + + +ops_h2 = [] +ops_h2.append([I(0)]) +ops_h2.append([I(0)]) +ops_h2.append( + [ + Y(0) @ Z(1) @ Y(2), + X(0) @ Z(1) @ X(2), + Y(1) @ Z(2) @ Y(3), + X(1) @ Z(2) @ X(3), + ] +) + +coeffs_h2_parity = [] +coeffs_h2_parity.append([0.0]) +coeffs_h2_parity.append([0.0]) +coeffs_h2_parity.append([-0.45445016, -0.45445016, -0.45445016, -0.45445016]) + + +ops_h2_parity = [] +ops_h2_parity.append([I(0)]) +ops_h2_parity.append([I(0)]) +ops_h2_parity.append( + [ + Y(0) @ Y(1), + X(0) @ X(1) @ Z(2), + Y(1) @ Y(2), + Z(0) @ X(1) @ X(2) @ Z(3), + ] +) + +h3p = ["H", "H", "H"] +x_h3p = np.array([[0.028, 0.054, 0.0], [0.986, 1.610, 0.0], [1.855, 0.002, 0.0]]) +coeffs_h3p = [] +coeffs_h3p.append( + [ + 0.47811232, + 0.47811232, + -0.39136385, + -0.39136385, + -0.39136385, + -0.39136385, + 0.26611147, + 0.26611147, + 0.26611147, + 0.26611147, + 0.71447791, + 0.71447791, + -0.11734959, + -0.11734959, + -0.11734959, + -0.11734959, + 0.24190978, + 0.24190978, + ] +) +coeffs_h3p.append( + [ + 0.27769368, + 0.27769368, + 0.26614699, + 0.26614699, + 0.26614699, + 0.26614699, + 0.39131162, + 0.39131162, + 0.39131162, + 0.39131162, + 0.16019825, + 0.16019825, + -0.23616713, + -0.23616713, + -0.23616713, + -0.23616713, + 0.39510807, + 0.39510807, + ] +) +coeffs_h3p.append([0.0]) + +ops_h3p = [] +ops_h3p.append( + [ + Z(0), + Z(1), + Y(0) @ Z(1) @ Y(2), + X(0) @ Z(1) @ X(2), + Y(1) @ Z(2) @ Y(3), + X(1) @ Z(2) @ X(3), + Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4), + X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4), + Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5), + X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5), + Z(2), + Z(3), + Y(2) @ Z(3) @ Y(4), + X(2) @ Z(3) @ X(4), + Y(3) @ Z(4) @ Y(5), + X(3) @ Z(4) @ X(5), + Z(4), + Z(5), + ] +) +ops_h3p.append( + [ + Z(0), + Z(1), + Y(0) @ Z(1) @ Y(2), + X(0) @ Z(1) @ X(2), + Y(1) @ Z(2) @ Y(3), + X(1) @ Z(2) @ X(3), + Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4), + X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4), + Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5), + X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5), + Z(2), + Z(3), + Y(2) @ Z(3) @ Y(4), + X(2) @ Z(3) @ X(4), + Y(3) @ Z(4) @ Y(5), + X(3) @ Z(4) @ X(5), + Z(4), + Z(5), + ] +) +ops_h3p.append([I(0)]) + + +h2o = ["H", "H", "O"] +x_h2o = np.array([0.0, 1.431, -0.887, 0.0, -1.431, -0.887, 0.0, 0.0, 0.222]) + +coeffs_h2o = [] +coeffs_h2o.append([-0.03700797, 0.03700797, 0.03700797, -0.03700797]) +coeffs_h2o.append([0.0]) +coeffs_h2o.append([0.28530461, 0.111, 0.111, -0.3710174, -0.3710174]) + +ops_h2o = [] +ops_h2o.append( + [ + X(0) @ Y(1) @ Y(2), + Y(0) @ Y(1) @ X(2), + Z(0) @ X(1) @ Z(3), + X(1) @ Z(2), + ] +) +ops_h2o.append([I(0)]) +ops_h2o.append( + [ + I(0), + Z(0), + Z(0) @ Z(1), + Z(2), + Z(1) @ Z(2) @ Z(3), + ] +) + + +@pytest.mark.parametrize( + ( + "symbols", + "geometry", + "charge", + "active_el", + "active_orb", + "mapping", + "coeffs", + "ops", + ), + [ + (h2, x_h2, 0, None, None, "jordan_wigner", coeffs_h2, ops_h2), + (h2, x_h2, 0, None, None, "parity", coeffs_h2_parity, ops_h2_parity), + (h3p, x_h3p, 1, None, None, "jordan_wigner", coeffs_h3p, ops_h3p), + (h2o, x_h2o, 0, 2, 2, "bravyi_kitaev", coeffs_h2o, ops_h2o), + ], +) +@pytest.mark.usefixtures("use_legacy_and_new_opmath") +def test_openfermion_molecular_dipole( + symbols, geometry, charge, active_el, active_orb, mapping, coeffs, ops, tol, tmpdir +): + r"""Test that molecular_dipole returns the correct dipole operator with openfermion backend.""" + + molecule = qml.qchem.Molecule(symbols, geometry, charge=charge) + dip = qml.qchem.molecular_dipole( + molecule, + method="openfermion", + active_electrons=active_el, + active_orbitals=active_orb, + mapping=mapping, + outpath=tmpdir.strpath, + ) + + assert len(dip) == len(ops) + + for i, _dip in enumerate(dip): + d_coeffs, d_ops = _dip.terms() + calc_coeffs = np.array(d_coeffs) + exp_coeffs = np.array(coeffs[i]) + assert np.allclose(calc_coeffs, exp_coeffs, **tol) + + r_ops = ops[i] + if not qml.operation.active_new_opmath(): + r_ops = [ + ( + qml.operation.Tensor(*obs.simplify()) + if isinstance(obs.simplify(), (qml.ops.op_math.Prod)) + else obs.simplify() + ) + for obs in ops[i] + ] + + assert all(isinstance(o1, o2.__class__) for o1, o2 in zip(d_ops, r_ops)) + assert all(qml.equal(o1, o2) for o1, o2 in zip(d_ops, r_ops)) + + +def test_molecular_dipole_error(): + r"""Test that molecular_dipole raises an error with unsupported backend and open-shell systems.""" + + 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', and 'openfermion' backends are supported"): + qchem.molecular_dipole(molecule, method="psi4") + + molecule = qchem.Molecule(symbols, geometry, mult=3) + with pytest.raises(ValueError, match="Open-shell systems are not supported"): + qchem.molecular_dipole(molecule) + + with pytest.raises(ValueError, match="'bksf' is not supported."): + qchem.molecular_dipole(molecule, mapping="bksf") + + +@pytest.mark.parametrize( + ("method", "args"), + [ + ( + "openfermion", + 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_dipole(method, args, tmpdir): + r"""Test that the generated operator 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) + + dipole = qchem.molecular_dipole( + molecule, + method=method, + args=args, + outpath=tmpdir.strpath, + ) + + assert all(np.isrealobj(op.terms()[0]) for op in dipole) + + +@pytest.mark.parametrize( + ("method"), + [ + "openfermion", + "dhf", + ], +) +def test_coordinate_units_for_molecular_dipole(method, tmpdir): + r"""Test that molecular_dipole generates correct operator 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") + dipole_bohr = qchem.molecular_dipole( + molecule_bohr, + method=method, + outpath=tmpdir.strpath, + ) + + molecule_ang = qchem.Molecule(symbols, geometry_ang, unit="angstrom") + dipole_ang = qchem.molecular_dipole( + molecule_ang, + method=method, + outpath=tmpdir.strpath, + ) + assert all(qml.equal(o1, o2) for o1, o2 in zip(dipole_ang, dipole_bohr)) diff --git a/tests/qchem/of_tests/test_observable_of.py b/tests/qchem/of_tests/test_observable_of.py index a74d37d4d0d..92adbf8796f 100644 --- a/tests/qchem/of_tests/test_observable_of.py +++ b/tests/qchem/of_tests/test_observable_of.py @@ -146,7 +146,7 @@ def test_observable(fermion_ops, init_term, mapping, terms_exp, custom_wires, mo msg1 = "Elements in the lists are expected to be of type 'FermionOperator'" -msg2 = "Please set 'mapping' to 'jordan_wigner' or 'bravyi_kitaev'" +msg2 = "Please set 'mapping' to 'jordan_wigner', 'parity', or 'bravyi_kitaev'" @pytest.mark.parametrize( From ffe2520617735a4215f031253978e7d207ef3414 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Thu, 30 May 2024 08:19:55 -0400 Subject: [PATCH 07/28] Added more tests --- doc/releases/changelog-dev.md | 3 + pennylane/qchem/openfermion_obs.py | 52 ++++++++++++++--- tests/qchem/of_tests/test_molecular_dipole.py | 56 +++++++++++++++++++ 3 files changed, 104 insertions(+), 7 deletions(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 0ac4a758940..dc062d8b6a6 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -18,6 +18,9 @@ * The sorting order of parameter-shift terms is now guaranteed to resolve ties in the absolute value with the sign of the shifts. [(#5582)](https://github.com/PennyLaneAI/pennylane/pull/5582) +* `molecular_dipole` function is added to calculate dipole operator using openfermion and dhf backends. + [(#5764)](https://github.com/PennyLaneAI/pennylane/pull/5764) +

Mid-circuit measurements and dynamic circuits

* The `dynamic_one_shot` transform uses a single auxiliary tape with a shot vector and `default.qubit` implements the loop over shots with `jax.vmap`. diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 0e41c6feda3..bed45550672 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -670,6 +670,48 @@ def molecular_dipole( cutoff=1.0e-16, ): # pylint:disable=too-many-arguments, too-many-statements r"""Generate the dipole moment operator for a molecule in the Pauli basis. + + The dipole operator in the second-quantized form is + + .. math:: + + \hat{D} = -\sum_{pq} d_{pq} [\hat{c}_{p\uparrow}^\dagger \hat{c}_{q\uparrow} + + \hat{c}_{p\downarrow}^\dagger \hat{c}_{q\downarrow}] - + \hat{D}_\mathrm{c} + \hat{D}_\mathrm{n}, + + where the matrix elements :math:`d_{pq}` are given by the integral of the position operator + :math:`\hat{{\bf r}}` over molecular orbitals :math:`\phi` + + .. math:: + + d_{pq} = \int \phi_p^*(r) \hat{{\bf r}} \phi_q(r) dr, + + and :math:`\hat{c}^{\dagger}` and :math:`\hat{c}` are the creation and annihilation operators, + respectively. The contribution of the core orbitals and nuclei are denoted by + :math:`\hat{D}_\mathrm{c}` and :math:`\hat{D}_\mathrm{n}`, respectively, which are computed as + + .. math:: + + \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii}, + + and + + .. math:: + + \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i, + + where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the + nuclear coordinates of the :math:`i`-th atom of the molecule. + + The fermonic dipole operator is then transformed to the qubit basis which gives + + .. math:: + + \hat{D} = \sum_{j} c_j P_j, + + where :math:`c_j` is a numerical coefficient and :math:`P_j` is a ternsor product of + single-qubit Pauli operators :math:`X, Y, Z, I`. + Args: molecule (~qchem.molecule.Molecule): the molecule object method (str): Quantum chemistry method used to solve the @@ -680,7 +722,7 @@ def molecular_dipole( are considered to be active. active_orbitals (int): Number of active orbitals. If not specified, all orbitals are considered to be active. - mapping (str): transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian + mapping (str): transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian. Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. outpath (str): path to the directory containing output files wires (Wires, list, tuple, dict): Custom wire mapping for connecting to Pennylane ansatz. For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label @@ -694,8 +736,8 @@ def molecular_dipole( Returns: list[pennylane.Hamiltonian]: the qubit observables corresponding to the components - :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator in - atomic units. + :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator. + """ @@ -712,15 +754,11 @@ def molecular_dipole( coordinates = molecule.coordinates if len(coordinates) == len(symbols) * 3: - print("1D coordinates") geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3)) geometry_hf = coordinates - print(geometry_dhf, geometry_hf) elif len(coordinates) == len(symbols): - print("2D coordinates") geometry_dhf = qml.numpy.array(coordinates) geometry_hf = coordinates.flatten() - print(geometry_dhf, geometry_hf) if molecule.mult != 1: raise ValueError( diff --git a/tests/qchem/of_tests/test_molecular_dipole.py b/tests/qchem/of_tests/test_molecular_dipole.py index 99cd4107b79..fe1309d76ab 100644 --- a/tests/qchem/of_tests/test_molecular_dipole.py +++ b/tests/qchem/of_tests/test_molecular_dipole.py @@ -61,6 +61,12 @@ ] ) +eig_h2 = [] +eig_h2.append([0., 0.]) +eig_h2.append([0., 0.]) +eig_h2.append([-1.81780062, -0.90890031, -0.90890031]) + + h3p = ["H", "H", "H"] x_h3p = np.array([[0.028, 0.054, 0.0], [0.986, 1.610, 0.0], [1.855, 0.002, 0.0]]) coeffs_h3p = [] @@ -157,6 +163,11 @@ ) ops_h3p.append([I(0)]) +eig_h3p = [] +eig_h3p.append([-3.15687028, -3.01293514, -3.01293514]) +eig_h3p.append([-2.00177188, -1.96904253, -1.96904253]) +eig_h3p.append([0., 0.]) + h2o = ["H", "H", "O"] x_h2o = np.array([0.0, 1.431, -0.887, 0.0, -1.431, -0.887, 0.0, 0.0, 0.222]) @@ -186,6 +197,10 @@ ] ) +eig_h2o=[] +eig_h2o.append([-0.14803188, -0.07401594, -0.07401594]) +eig_h2o.append([0., 0.]) +eig_h2o.append([-0.67873019, -0.45673019, -0.45673019]) @pytest.mark.parametrize( ( @@ -243,6 +258,47 @@ def test_openfermion_molecular_dipole( assert all(isinstance(o1, o2.__class__) for o1, o2 in zip(d_ops, r_ops)) assert all(qml.equal(o1, o2) for o1, o2 in zip(d_ops, r_ops)) +@pytest.mark.parametrize( + ( + "symbols", + "geometry", + "charge", + "active_el", + "active_orb", + "mapping", + "eig_ref", + ), + [ + (h2, x_h2, 0, None, None, "jordan_wigner", eig_h2), + (h2, x_h2, 0, None, None, "parity", eig_h2), + (h3p, x_h3p, 1, None, None, "jordan_wigner", eig_h3p), + (h2o, x_h2o, 0, 2, 2, "bravyi_kitaev", eig_h2o), + ], +) +@pytest.mark.usefixtures("use_legacy_and_new_opmath") +def test_differentiable_molecular_dipole( + symbols, geometry, charge, active_el, active_orb, mapping, eig_ref, tol, tmpdir +): + r"""Test that molecular_dipole returns the correct eigenvalues with the dhf backend.""" + + molecule = qml.qchem.Molecule(symbols, geometry, charge=charge) + dip_dhf = qml.qchem.molecular_dipole( + molecule, + method="dhf", + active_electrons=active_el, + active_orbitals=active_orb, + mapping=mapping, + outpath=tmpdir.strpath, + ) + + for idx, dip in enumerate(dip_dhf): + wires = dip.wires + if not wires: + eig =[0, 0] + else: + eig = qml.eigvals(qml.SparseHamiltonian(dip.sparse_matrix(), wires = wires), k=3) + assert(np.allclose(np.sort(eig), np.sort(eig_ref[idx]))) + def test_molecular_dipole_error(): r"""Test that molecular_dipole raises an error with unsupported backend and open-shell systems.""" From 062c2a6900946f3d5184419a888f08a5f5f881b5 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Thu, 30 May 2024 08:49:59 -0400 Subject: [PATCH 08/28] Fixed CI --- tests/qchem/of_tests/test_molecular_dipole.py | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/tests/qchem/of_tests/test_molecular_dipole.py b/tests/qchem/of_tests/test_molecular_dipole.py index fe1309d76ab..ce81e3d0ca0 100644 --- a/tests/qchem/of_tests/test_molecular_dipole.py +++ b/tests/qchem/of_tests/test_molecular_dipole.py @@ -21,7 +21,6 @@ from pennylane import I, X, Y, Z from pennylane import numpy as np from pennylane import qchem -from pennylane.operation import active_new_opmath h2 = ["H", "H"] x_h2 = np.array([0.0, 0.0, -0.661, 0.0, 0.0, 0.661]) @@ -62,8 +61,8 @@ ) eig_h2 = [] -eig_h2.append([0., 0.]) -eig_h2.append([0., 0.]) +eig_h2.append([0.0, 0.0]) +eig_h2.append([0.0, 0.0]) eig_h2.append([-1.81780062, -0.90890031, -0.90890031]) @@ -166,7 +165,7 @@ eig_h3p = [] eig_h3p.append([-3.15687028, -3.01293514, -3.01293514]) eig_h3p.append([-2.00177188, -1.96904253, -1.96904253]) -eig_h3p.append([0., 0.]) +eig_h3p.append([0.0, 0.0]) h2o = ["H", "H", "O"] @@ -197,11 +196,12 @@ ] ) -eig_h2o=[] +eig_h2o = [] eig_h2o.append([-0.14803188, -0.07401594, -0.07401594]) -eig_h2o.append([0., 0.]) +eig_h2o.append([0.0, 0.0]) eig_h2o.append([-0.67873019, -0.45673019, -0.45673019]) + @pytest.mark.parametrize( ( "symbols", @@ -258,6 +258,7 @@ def test_openfermion_molecular_dipole( assert all(isinstance(o1, o2.__class__) for o1, o2 in zip(d_ops, r_ops)) assert all(qml.equal(o1, o2) for o1, o2 in zip(d_ops, r_ops)) + @pytest.mark.parametrize( ( "symbols", @@ -277,7 +278,7 @@ def test_openfermion_molecular_dipole( ) @pytest.mark.usefixtures("use_legacy_and_new_opmath") def test_differentiable_molecular_dipole( - symbols, geometry, charge, active_el, active_orb, mapping, eig_ref, tol, tmpdir + symbols, geometry, charge, active_el, active_orb, mapping, eig_ref, tmpdir ): r"""Test that molecular_dipole returns the correct eigenvalues with the dhf backend.""" @@ -294,11 +295,11 @@ def test_differentiable_molecular_dipole( for idx, dip in enumerate(dip_dhf): wires = dip.wires if not wires: - eig =[0, 0] + eig = [0, 0] else: - eig = qml.eigvals(qml.SparseHamiltonian(dip.sparse_matrix(), wires = wires), k=3) - assert(np.allclose(np.sort(eig), np.sort(eig_ref[idx]))) - + eig = qml.eigvals(qml.SparseHamiltonian(dip.sparse_matrix(), wires=wires), k=3) + assert np.allclose(np.sort(eig), np.sort(eig_ref[idx])) + def test_molecular_dipole_error(): r"""Test that molecular_dipole raises an error with unsupported backend and open-shell systems.""" From 65e037a20ec8950bf1c78cc3e4adaeed994509c2 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Thu, 30 May 2024 09:34:24 -0400 Subject: [PATCH 09/28] [skip ci] resolved conflicts --- doc/code/qml_qchem.rst | 3 +-- doc/releases/changelog-dev.md | 10 +++------- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/doc/code/qml_qchem.rst b/doc/code/qml_qchem.rst index 68147b590a9..c5e86986e90 100644 --- a/doc/code/qml_qchem.rst +++ b/doc/code/qml_qchem.rst @@ -185,7 +185,6 @@ in ``molecular_hamiltonian``: symbols = ["H", "H"] geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]]) - molecule = qml.qchem.Molecule(symbols, geometry, charge=0, mult=1, basis_name='sto-3g') hamiltonian, qubits = qml.qchem.molecular_hamiltonian(molecule, method='pyscf') @@ -195,4 +194,4 @@ the user with .. code-block:: bash pip install openfermionpyscf - pip install pyscf \ No newline at end of file + pip install pyscf diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 820e9a45bb6..cb7d485f256 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -21,9 +21,6 @@ * The sorting order of parameter-shift terms is now guaranteed to resolve ties in the absolute value with the sign of the shifts. [(#5582)](https://github.com/PennyLaneAI/pennylane/pull/5582) -* `molecular_dipole` function is added to calculate dipole operator using openfermion and dhf backends. - [(#5764)](https://github.com/PennyLaneAI/pennylane/pull/5764) -

Mid-circuit measurements and dynamic circuits

* The `dynamic_one_shot` transform uses a single auxiliary tape with a shot vector and `default.qubit` implements the loop over shots with `jax.vmap`. @@ -86,10 +83,6 @@ `qml.devices.Device`, which follows the new device API. [(#5581)](https://github.com/PennyLaneAI/pennylane/pull/5581) -* `qml.qchem.Molecule` and qchem functions that take `Molecule` object as an argument now work with coordinates - provided in units Angstrom along with Bohr. - [(#5694)](https://github.com/PennyLaneAI/pennylane/pull/5694) - * The `dtype` for `eigvals` of `X`, `Y`, `Z` and `Hadamard` is changed from `int` to `float`, making them consistent with the other observables. The `dtype` of the returned values when sampling these observables (e.g. `qml.sample(X(0))`) is also changed to `float`. @@ -137,6 +130,9 @@ [(#5758)](https://github.com/PennyLaneAI/pennylane/pull/5758/) [(#5638)](https://github.com/PennyLaneAI/pennylane/pull/5638/) +* `molecular_dipole` function is added to calculate dipole operator using openfermion and dhf backends. + [(#5764)](https://github.com/PennyLaneAI/pennylane/pull/5764) +

Community contributions 🥳

* Implemented kwargs (`check_interface`, `check_trainability`, `rtol` and `atol`) support in `qml.equal` for the operators `Pow`, `Adjoint`, `Exp`, and `SProd`. From 05f834631720b4db5baff89604cefd1d431b7822 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 09:57:05 -0400 Subject: [PATCH 10/28] Update doc/releases/changelog-dev.md Co-authored-by: Utkarsh --- doc/releases/changelog-dev.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index cb7d485f256..6bfd33b4005 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -130,7 +130,7 @@ [(#5758)](https://github.com/PennyLaneAI/pennylane/pull/5758/) [(#5638)](https://github.com/PennyLaneAI/pennylane/pull/5638/) -* `molecular_dipole` function is added to calculate dipole operator using openfermion and dhf backends. +* `qml.qchem.molecular_dipole` function is added for calculating the dipole operator using "dhf" and "openfermion" backends. [(#5764)](https://github.com/PennyLaneAI/pennylane/pull/5764)

Community contributions 🥳

From 6528379698d533bef4c25465919c6555e460950e Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 09:57:41 -0400 Subject: [PATCH 11/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 01b9677a639..388a72c1603 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -691,13 +691,7 @@ def molecular_dipole( :math:`\hat{D}_\mathrm{c}` and :math:`\hat{D}_\mathrm{n}`, respectively, which are computed as .. math:: - - \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii}, - - and - - .. math:: - + \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii} \quad \text{and} \quad \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i, where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the From 4676911c55fcf61a7e5e44325cc25924a70cfb45 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 09:57:47 -0400 Subject: [PATCH 12/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 388a72c1603..de965a4b398 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -703,7 +703,7 @@ def molecular_dipole( \hat{D} = \sum_{j} c_j P_j, - where :math:`c_j` is a numerical coefficient and :math:`P_j` is a ternsor product of + where :math:`c_j` is a numerical coefficient and :math:`P_j` is a tensor product of single-qubit Pauli operators :math:`X, Y, Z, I`. Args: From 1a502bc30b74c695c13098d332cbee53d4808567 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 09:57:56 -0400 Subject: [PATCH 13/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index de965a4b398..808bf9c5562 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -697,7 +697,7 @@ def molecular_dipole( where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the nuclear coordinates of the :math:`i`-th atom of the molecule. - The fermonic dipole operator is then transformed to the qubit basis which gives + The fermionic dipole operator is then transformed to the qubit basis, which gives .. math:: From be368309ec218ff1626ddf015b4de492244d61d4 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 09:58:03 -0400 Subject: [PATCH 14/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 808bf9c5562..416f51192d9 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -687,7 +687,7 @@ def molecular_dipole( d_{pq} = \int \phi_p^*(r) \hat{{\bf r}} \phi_q(r) dr, and :math:`\hat{c}^{\dagger}` and :math:`\hat{c}` are the creation and annihilation operators, - respectively. The contribution of the core orbitals and nuclei are denoted by + respectively. The contribution of the core orbitals and nuclei is denoted by :math:`\hat{D}_\mathrm{c}` and :math:`\hat{D}_\mathrm{n}`, respectively, which are computed as .. math:: From c534e647ffd572d8d5288ecf34a400d997177c3c Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 10:04:19 -0400 Subject: [PATCH 15/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 416f51192d9..a8a6b066a27 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -787,7 +787,7 @@ def molecular_dipole( requires_grad = args is not None dip = ( qml.qchem.dipole_moment(mol, cutoff=cutoff, core=core, active=active, mapping=mapping)( - *args + *(args or []) ) if requires_grad else qml.qchem.dipole_moment( From 4e7bb867fc3cf3f7bf22a5391fa797398f3eff02 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 10:11:19 -0400 Subject: [PATCH 16/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index a8a6b066a27..09d48021eca 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -710,7 +710,7 @@ def molecular_dipole( 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, `or to + to specify the built-in differentiable Hartree-Fock solver, 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. From f7b88807bfa580e908aec773d4e58d0d0822db97 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Mon, 3 Jun 2024 10:14:15 -0400 Subject: [PATCH 17/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 09d48021eca..349e4215374 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -747,7 +747,7 @@ def molecular_dipole( symbols = molecule.symbols coordinates = molecule.coordinates - if len(coordinates) == len(symbols) * 3: + if qml.math.shape(coordinates)[0] == len(symbols) * 3: geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3)) geometry_hf = coordinates elif len(coordinates) == len(symbols): From e7713164cfd6309c421996608110590e568adde3 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Tue, 4 Jun 2024 09:25:40 -0400 Subject: [PATCH 18/28] [skip ci] Addressed review comments --- pennylane/qchem/openfermion_obs.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 349e4215374..30104c83216 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -687,11 +687,11 @@ def molecular_dipole( d_{pq} = \int \phi_p^*(r) \hat{{\bf r}} \phi_q(r) dr, and :math:`\hat{c}^{\dagger}` and :math:`\hat{c}` are the creation and annihilation operators, - respectively. The contribution of the core orbitals and nuclei is denoted by + respectively. The contribution of the core orbitals and nuclei are denoted by :math:`\hat{D}_\mathrm{c}` and :math:`\hat{D}_\mathrm{n}`, respectively, which are computed as .. math:: - \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii} \quad \text{and} \quad + \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii} \quad \text{and} \quad \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i, where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the @@ -704,10 +704,12 @@ def molecular_dipole( \hat{D} = \sum_{j} c_j P_j, where :math:`c_j` is a numerical coefficient and :math:`P_j` is a tensor product of - single-qubit Pauli operators :math:`X, Y, Z, I`. + single-qubit Pauli operators :math:`X, Y, Z, I`. The qubit observables corresponding + to the components :math:`\hat{D}_x`, :math:`\hat{D}_y`, and :math:`\hat{D}_z` of the + dipole operator are then computed separately. Args: - molecule (~qchem.molecule.Molecule): the molecule object + 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, or ``method="openfermion"`` to @@ -716,20 +718,22 @@ def molecular_dipole( are considered to be active. active_orbitals (int): Number of active orbitals. If not specified, all orbitals are considered to be active. - mapping (str): transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian. Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. - outpath (str): path to the directory containing output files - wires (Wires, list, tuple, dict): Custom wire mapping for connecting to Pennylane ansatz. + mapping (str): Transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian. + Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. + outpath (str): Path to the directory containing output files + wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator to + an observable measurable in a Pennylane ansatz. For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label 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. - args (array[array[float]]): initial values of the differentiable parameters + args (array[array[float]]): Initial values of the differentiable parameters cutoff (float): Cutoff value for including the matrix elements :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements with absolute value less than ``cutoff`` are neglected. Returns: - list[pennylane.Hamiltonian]: the qubit observables corresponding to the components + list[pennylane.Hamiltonian]: The qubit observables corresponding to the components :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator. From a9d78a1f9d834ea72cdf83a6e8ff4ff3e21179ba Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Tue, 4 Jun 2024 09:31:04 -0400 Subject: [PATCH 19/28] [skip ci] Added example --- pennylane/qchem/openfermion_obs.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 30104c83216..f9370048f75 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -736,7 +736,33 @@ def molecular_dipole( list[pennylane.Hamiltonian]: The qubit observables corresponding to the components :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator. + **Example** + >>> symbols = ["H", "H", "H"] + >>> coordinates = np.array([[0.028, 0.054, 0.0], [0.986, 1.610, 0.0], [1.855, 0.002, 0.0]]) + >>> mol = qml.qchem.Molecule(symbols, coordinates, charge=1) + >>> dipole_obs = qml.qchem.molecular_dipole(mol, method="openfermion") + >>> dipole_obs[0] # x-component of D + ( + 0.4781123173263876 * Z(0) + + 0.4781123173263876 * Z(1) + + -0.3913638489489803 * (Y(0) @ Z(1) @ Y(2)) + + -0.3913638489489803 * (X(0) @ Z(1) @ X(2)) + + -0.3913638489489803 * (Y(1) @ Z(2) @ Y(3)) + + -0.3913638489489803 * (X(1) @ Z(2) @ X(3)) + + 0.2661114704527088 * (Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4)) + + 0.2661114704527088 * (X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4)) + + 0.2661114704527088 * (Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5)) + + 0.2661114704527088 * (X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5)) + + 0.7144779061810713 * Z(2) + + 0.7144779061810713 * Z(3) + + -0.11734958781031017 * (Y(2) @ Z(3) @ Y(4)) + + -0.11734958781031017 * (X(2) @ Z(3) @ X(4)) + + -0.11734958781031017 * (Y(3) @ Z(4) @ Y(5)) + + -0.11734958781031017 * (X(3) @ Z(4) @ X(5)) + + 0.24190977644645698 * Z(4) + + 0.24190977644645698 * Z(5) + ) """ if method not in ["dhf", "openfermion"]: From 480298b68c91dcd4a71a80f33e38e9edb539c51f Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Wed, 5 Jun 2024 08:16:07 -0400 Subject: [PATCH 20/28] Increased test coverage --- pennylane/qchem/openfermion_obs.py | 3 ++ tests/qchem/of_tests/test_molecular_dipole.py | 34 +++++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index f9370048f75..3dea709db8c 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -824,6 +824,9 @@ def molecular_dipole( mol, cutoff=cutoff, core=core, active=active, mapping=mapping )() ) + if wires: + dip = [qml.map_wires(op, wires_map) for op in dip] + return dip dip = qml.qchem.dipole_of( diff --git a/tests/qchem/of_tests/test_molecular_dipole.py b/tests/qchem/of_tests/test_molecular_dipole.py index ce81e3d0ca0..52b96023c6a 100644 --- a/tests/qchem/of_tests/test_molecular_dipole.py +++ b/tests/qchem/of_tests/test_molecular_dipole.py @@ -300,6 +300,40 @@ def test_differentiable_molecular_dipole( eig = qml.eigvals(qml.SparseHamiltonian(dip.sparse_matrix(), wires=wires), k=3) assert np.allclose(np.sort(eig), np.sort(eig_ref[idx])) +@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_dipole(wiremap, tmpdir): + r"""Test that the generated dipole operator 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 = "dhf" + + mol = qchem.Molecule(symbols, geometry) + dipole_dhf = qchem.molecular_dipole( + mol, + method=method, + wires=wiremap, + outpath=tmpdir.strpath, + ) + + method = "openfermion" + dipole_of = qchem.molecular_dipole( + mol, + method=method, + wires=wiremap, + outpath=tmpdir.strpath, + ) + + assert all(set(dip.wires).issubset(set(wiremap)) for dip in dipole_dhf) + assert all(set(dip.wires).issubset(set(wiremap)) for dip in dipole_of) def test_molecular_dipole_error(): r"""Test that molecular_dipole raises an error with unsupported backend and open-shell systems.""" From 8d8679f322a5009d286ccad728e46b67f77c708d Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Wed, 5 Jun 2024 10:28:47 -0400 Subject: [PATCH 21/28] Fixed CI --- tests/qchem/of_tests/test_molecular_dipole.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/qchem/of_tests/test_molecular_dipole.py b/tests/qchem/of_tests/test_molecular_dipole.py index 52b96023c6a..29aad2b2222 100644 --- a/tests/qchem/of_tests/test_molecular_dipole.py +++ b/tests/qchem/of_tests/test_molecular_dipole.py @@ -300,6 +300,7 @@ def test_differentiable_molecular_dipole( eig = qml.eigvals(qml.SparseHamiltonian(dip.sparse_matrix(), wires=wires), k=3) assert np.allclose(np.sort(eig), np.sort(eig_ref[idx])) + @pytest.mark.usefixtures("use_legacy_and_new_opmath") @pytest.mark.parametrize( ("wiremap"), @@ -335,6 +336,7 @@ def test_custom_wiremap_dipole(wiremap, tmpdir): assert all(set(dip.wires).issubset(set(wiremap)) for dip in dipole_dhf) assert all(set(dip.wires).issubset(set(wiremap)) for dip in dipole_of) + def test_molecular_dipole_error(): r"""Test that molecular_dipole raises an error with unsupported backend and open-shell systems.""" From e47aa6425608b87df0144cf27777acfa44d6cd7c Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Wed, 5 Jun 2024 12:46:30 -0400 Subject: [PATCH 22/28] Addressed reviews --- pennylane/qchem/dipole.py | 2 +- pennylane/qchem/openfermion_obs.py | 11 +++-------- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/pennylane/qchem/dipole.py b/pennylane/qchem/dipole.py index d4c36a7505c..05c04a35ef0 100644 --- a/pennylane/qchem/dipole.py +++ b/pennylane/qchem/dipole.py @@ -277,7 +277,7 @@ def dipole_moment(mol, cutoff=1.0e-16, core=None, active=None, mapping="jordan_w cutoff (float): cutoff value for discarding the negligible dipole moment integrals core (list[int]): indices of the core orbitals active (list[int]): indices of the active orbitals - mapping (str): Specifies the transformation to map the fermionic Hamiltonian to the + mapping (str): Specifies the transformation to map the fermionic dipole operator to the Pauli basis. Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. Returns: diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index 3dea709db8c..cd6a3cac669 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -789,12 +789,6 @@ def molecular_dipole( "Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule." ) - 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))) - core, active = qml.qchem.active_space( molecule.n_electrons, molecule.n_orbitals, molecule.mult, active_electrons, active_orbitals ) @@ -824,7 +818,10 @@ def molecular_dipole( mol, cutoff=cutoff, core=core, active=active, mapping=mapping )() ) + 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))) dip = [qml.map_wires(op, wires_map) for op in dip] return dip @@ -1001,7 +998,6 @@ def decompose(hf_file, mapping="jordan_wigner", core=None, active=None): 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 @@ -1057,7 +1053,6 @@ def molecular_hamiltonian(*args, **kwargs): 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 = ['H', 'H'] From 1c0c19504c78109ecfb6a4d788c2f8b93b1e678a Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Wed, 5 Jun 2024 12:52:37 -0400 Subject: [PATCH 23/28] Addressed reviews --- tests/qchem/of_tests/test_molecular_dipole.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/qchem/of_tests/test_molecular_dipole.py b/tests/qchem/of_tests/test_molecular_dipole.py index 29aad2b2222..070057b29ad 100644 --- a/tests/qchem/of_tests/test_molecular_dipole.py +++ b/tests/qchem/of_tests/test_molecular_dipole.py @@ -24,6 +24,7 @@ h2 = ["H", "H"] x_h2 = np.array([0.0, 0.0, -0.661, 0.0, 0.0, 0.661]) +# Reference dipole operator coeffs_h2 = [] coeffs_h2.append([0.0]) coeffs_h2.append([0.0]) @@ -42,6 +43,7 @@ ] ) +# Reference dipole operator coeffs_h2_parity = [] coeffs_h2_parity.append([0.0]) coeffs_h2_parity.append([0.0]) @@ -60,6 +62,7 @@ ] ) +# Reference eigenvalues eig_h2 = [] eig_h2.append([0.0, 0.0]) eig_h2.append([0.0, 0.0]) @@ -68,6 +71,7 @@ h3p = ["H", "H", "H"] x_h3p = np.array([[0.028, 0.054, 0.0], [0.986, 1.610, 0.0], [1.855, 0.002, 0.0]]) +# Reference dipole operator coeffs_h3p = [] coeffs_h3p.append( [ @@ -162,6 +166,7 @@ ) ops_h3p.append([I(0)]) +# Reference eigenvalues eig_h3p = [] eig_h3p.append([-3.15687028, -3.01293514, -3.01293514]) eig_h3p.append([-2.00177188, -1.96904253, -1.96904253]) @@ -171,6 +176,7 @@ h2o = ["H", "H", "O"] x_h2o = np.array([0.0, 1.431, -0.887, 0.0, -1.431, -0.887, 0.0, 0.0, 0.222]) +# Reference dipole operator coeffs_h2o = [] coeffs_h2o.append([-0.03700797, 0.03700797, 0.03700797, -0.03700797]) coeffs_h2o.append([0.0]) @@ -196,6 +202,7 @@ ] ) +# Reference eigenvalues eig_h2o = [] eig_h2o.append([-0.14803188, -0.07401594, -0.07401594]) eig_h2o.append([0.0, 0.0]) From b6b094449a490be1e0842bc631dfcce6c3914729 Mon Sep 17 00:00:00 2001 From: Diksha Dhawan <40900030+ddhawan11@users.noreply.github.com> Date: Thu, 6 Jun 2024 10:02:12 -0400 Subject: [PATCH 24/28] Update pennylane/qchem/openfermion_obs.py Co-authored-by: Utkarsh --- pennylane/qchem/openfermion_obs.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py index cd6a3cac669..70a2c078d63 100644 --- a/pennylane/qchem/openfermion_obs.py +++ b/pennylane/qchem/openfermion_obs.py @@ -818,7 +818,6 @@ def molecular_dipole( mol, cutoff=cutoff, core=core, active=active, mapping=mapping )() ) - 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))) From 9b0b5fb959473b988f4c9cbf3a63c28d5d0ba463 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Thu, 6 Jun 2024 12:23:41 -0400 Subject: [PATCH 25/28] Resolved conflicts --- pennylane/qchem/openfermion_pyscf.py | 1091 ++++++++++++++++++++++++++ 1 file changed, 1091 insertions(+) create mode 100644 pennylane/qchem/openfermion_pyscf.py diff --git a/pennylane/qchem/openfermion_pyscf.py b/pennylane/qchem/openfermion_pyscf.py new file mode 100644 index 00000000000..bc396b37fbe --- /dev/null +++ b/pennylane/qchem/openfermion_pyscf.py @@ -0,0 +1,1091 @@ +# Copyright 2018-2020 Xanadu Quantum Technologies Inc. + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""This module contains functions to construct many-body observables with ``OpenFermion-PySCF``. +""" +# pylint: disable=too-many-arguments, too-few-public-methods, too-many-branches, unused-variable +# pylint: disable=consider-using-generator, protected-access +import os + +import numpy as np + +import pennylane as qml + +from .basis_data import atomic_numbers + +# Bohr-Angstrom correlation coefficient (https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0) +bohr_angs = 0.529177210903 + + +def _import_of(): + """Import openfermion and openfermionpyscf.""" + try: + # pylint: disable=import-outside-toplevel, unused-import, multiple-imports + import openfermion + import openfermionpyscf + except ImportError as Error: + raise ImportError( + "This feature requires openfermionpyscf. " + "It can be installed with: pip install openfermionpyscf." + ) from Error + + 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. + + The second-quantized operator of the fermionic many-body system can combine one-particle + and two-particle operators as in the case of electronic Hamiltonians :math:`\hat{H}`: + + .. math:: + + \hat{H} = \sum_{\alpha, \beta} \langle \alpha \vert \hat{t}^{(1)} + + \cdots + \hat{t}^{(n)} \vert \beta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta + + \frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta} + \langle \alpha, \beta \vert \hat{v}^{(1)} + \cdots + \hat{v}^{(n)} + \vert \gamma, \delta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta^\dagger + \hat{c}_\gamma \hat{c}_\delta + + In the latter equations the indices :math:`\alpha, \beta, \gamma, \delta` run over the + basis of single-particle states. The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` + are the particle creation and annihilation operators, respectively. + :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` denotes the matrix element of + the single-particle operator :math:`\hat{t}` entering the observable. For example, + in electronic structure calculations, this is the case for the kinetic energy operator, + the nuclei Coulomb potential, or any other external fields included in the Hamiltonian. + On the other hand, :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` + denotes the matrix element of the two-particle operator :math:`\hat{v}`, for example, the + Coulomb interaction between the electrons. + + - The observable is built by adding the operators + :math:`\sum_{\alpha, \beta} t_{\alpha\beta}^{(i)} + \hat{c}_\alpha^\dagger \hat{c}_\beta` and + :math:`\frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta} + v_{\alpha\beta\gamma\delta}^{(i)} + \hat{c}_\alpha^\dagger \hat{c}_\beta^\dagger \hat{c}_\gamma \hat{c}_\delta`. + + - Second-quantized operators contributing to the + many-body observable must be represented using the `FermionOperator + `_ data structure as implemented in OpenFermion. + See the functions :func:`~.one_particle` and :func:`~.two_particle` to build the + FermionOperator representations of one-particle and two-particle operators. + + - The function uses tools of `OpenFermion `_ + to map the resulting fermionic Hamiltonian to the basis of Pauli matrices via the + Jordan-Wigner or Bravyi-Kitaev transformation. Finally, the qubit operator is converted + to a PennyLane observable by the function :func:`~.convert_observable`. + + Args: + fermion_ops (list[FermionOperator]): list containing the FermionOperator data structures + representing the one-particle and/or two-particle operators entering the many-body + observable + init_term (float): Any quantity required to initialize the many-body observable. For + example, this can be used to pass the nuclear-nuclear repulsion energy :math:`V_{nn}` + which is typically included in the electronic Hamiltonian of molecules. + mapping (str): Specifies the fermion-to-qubit mapping. Input values can + be ``'jordan_wigner'``, ``'parity'``, or ``'bravyi_kitaev'``. + wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator + to an observable measurable in a PennyLane ansatz. + For types Wires/list/tuple, each item in the iterable represents a wire label + corresponding to the qubit number equal to its index. + For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted. + If None, will use identity map (e.g. 0->0, 1->1, ...). + + Returns: + pennylane.Hamiltonian: the fermionic-to-qubit transformed observable + + **Example** + + >>> t = FermionOperator("0^ 0", 0.5) + FermionOperator("1^ 1", 0.25) + >>> v = FermionOperator("1^ 0^ 0 1", -0.15) + FermionOperator("2^ 0^ 2 0", 0.3) + >>> observable([t, v], mapping="jordan_wigner") + ( + 0.2625 * I(0) + + -0.1375 * Z(0) + + -0.0875 * Z(1) + + -0.0375 * (Z(0) @ Z(1)) + + 0.075 * Z(2) + + -0.075 * (Z(0) @ Z(2)) + ) + """ + openfermion, _ = _import_of() + + if mapping.strip().lower() not in ("jordan_wigner", "parity", "bravyi_kitaev"): + raise TypeError( + f"The '{mapping}' transformation is not available. \n " + f"Please set 'mapping' to 'jordan_wigner', 'parity', or 'bravyi_kitaev'." + ) + + # Initialize the FermionOperator + mb_obs = openfermion.ops.FermionOperator("") * init_term + for ops in fermion_ops: + if not isinstance(ops, openfermion.ops.FermionOperator): + raise TypeError( + f"Elements in the lists are expected to be of type 'FermionOperator'; got {type(ops)}" + ) + mb_obs += ops + + # Map the fermionic operator to a qubit operator + if mapping.strip().lower() == "bravyi_kitaev": + return qml.qchem.convert.import_operator( + openfermion.transforms.bravyi_kitaev(mb_obs), wires=wires + ) + + if mapping.strip().lower() == "parity": + qubits = openfermion.count_qubits(mb_obs) + if qubits == 0: + return 0.0 * qml.I(0) + binary_code = openfermion.parity_code(qubits) + return qml.qchem.convert.import_operator( + openfermion.transforms.binary_code_transform(mb_obs, binary_code), wires=wires + ) + + return qml.qchem.convert.import_operator( + openfermion.transforms.jordan_wigner(mb_obs), wires=wires + ) + + +def one_particle(matrix_elements, core=None, active=None, cutoff=1.0e-12): + r"""Generates the `FermionOperator `_ representing a given one-particle operator + required to build many-body qubit observables. + + Second quantized one-particle operators are expanded in the basis of single-particle + states as + + .. math:: + + \hat{T} = \sum_{\alpha, \beta} \langle \alpha \vert \hat{t} \vert \beta \rangle + [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} + + \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}]. + + In the equation above the indices :math:`\alpha, \beta` run over the basis of spatial + orbitals :math:`\phi_\alpha(r)`. Since the operator :math:`\hat{t}` acts only on the + spatial coordinates, the spin quantum numbers are indicated explicitly with the up/down arrows. + The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation + and annihilation operators, respectively, and + :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` denotes the matrix elements of + the operator :math:`\hat{t}` + + .. math:: + + \langle \alpha \vert \hat{t} \vert \beta \rangle = \int dr ~ \phi_\alpha^*(r) + \hat{t}(r) \phi_\beta(r). + + If an active space is defined (see :func:`~.active_space`), the summation indices + run over the active orbitals and the contribution due to core orbitals is computed as + :math:`t_\mathrm{core} = 2 \sum_{\alpha\in \mathrm{core}} + \langle \alpha \vert \hat{t} \vert \alpha \rangle`. + + Args: + matrix_elements (array[float]): 2D NumPy array with the matrix elements + :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` + core (list): indices of core orbitals, i.e., the orbitals that are + not correlated in the many-body wave function + active (list): indices of active orbitals, i.e., the orbitals used to + build the correlated many-body wave function + cutoff (float): Cutoff value for including matrix elements. The + matrix elements with absolute value less than ``cutoff`` are neglected. + + Returns: + FermionOperator: an instance of OpenFermion's ``FermionOperator`` representing the + one-particle operator :math:`\hat{T}`. + + **Example** + + >>> import numpy as np + >>> matrix_elements = np.array([[-1.27785301e+00, 0.00000000e+00], + ... [ 1.52655666e-16, -4.48299696e-01]]) + >>> t_op = one_particle(matrix_elements) + >>> print(t_op) + -1.277853006156875 [0^ 0] + + -1.277853006156875 [1^ 1] + + -0.44829969610163756 [2^ 2] + + -0.44829969610163756 [3^ 3] + """ + openfermion, _ = _import_of() + + orbitals = matrix_elements.shape[0] + + if matrix_elements.ndim != 2: + raise ValueError( + f"'matrix_elements' must be a 2D array; got matrix_elements.ndim = {matrix_elements.ndim}" + ) + + if not core: + t_core = 0 + else: + if any(i > orbitals - 1 or i < 0 for i in core): + raise ValueError( + f"Indices of core orbitals must be between 0 and {orbitals - 1}; got core = {core}" + ) + + # Compute contribution due to core orbitals + t_core = 2 * np.sum(matrix_elements[np.array(core), np.array(core)]) + + if active is None: + if not core: + active = list(range(orbitals)) + else: + active = [i for i in range(orbitals) if i not in core] + + if any(i > orbitals - 1 or i < 0 for i in active): + raise ValueError( + f"Indices of active orbitals must be between 0 and {orbitals - 1}; got active = {active}" + ) + + # Indices of the matrix elements with absolute values >= cutoff + indices = np.nonzero(np.abs(matrix_elements) >= cutoff) + + # Single out the indices of active orbitals + num_indices = len(indices[0]) + pairs = [ + [indices[0][i], indices[1][i]] + for i in range(num_indices) + if all(indices[j][i] in active for j in range(len(indices))) + ] + + # Build the FermionOperator representing T + t_op = openfermion.ops.FermionOperator("") * t_core + for pair in pairs: + alpha, beta = pair + element = matrix_elements[alpha, beta] + + # spin-up term + a = 2 * active.index(alpha) + b = 2 * active.index(beta) + t_op += openfermion.ops.FermionOperator(((a, 1), (b, 0)), element) + + # spin-down term + t_op += openfermion.ops.FermionOperator(((a + 1, 1), (b + 1, 0)), element) + + return t_op + + +def two_particle(matrix_elements, core=None, active=None, cutoff=1.0e-12): + r"""Generates the `FermionOperator `_ representing a given two-particle operator + required to build many-body qubit observables. + + Second quantized two-particle operators are expanded in the basis of single-particle + states as + + .. math:: + + \hat{V} = \frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta} + \langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle + ~ &[& \hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow}^\dagger + \hat{c}_{\gamma\uparrow} \hat{c}_{\delta\uparrow} + \hat{c}_{\alpha\uparrow}^\dagger + \hat{c}_{\beta\downarrow}^\dagger \hat{c}_{\gamma\downarrow} \hat{c}_{\delta\uparrow} \\ + &+& \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\uparrow}^\dagger + \hat{c}_{\gamma\uparrow} \hat{c}_{\delta\downarrow} + \hat{c}_{\alpha\downarrow}^\dagger + \hat{c}_{\beta\downarrow}^\dagger \hat{c}_{\gamma\downarrow} \hat{c}_{\delta\downarrow}~]. + + In the equation above the indices :math:`\alpha, \beta, \gamma, \delta` run over the basis + of spatial orbitals :math:`\phi_\alpha(r)`. Since the operator :math:`v` acts only on the + spatial coordinates the spin quantum numbers are indicated explicitly with the up/down arrows. + The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation and + annihilation operators, respectively, and + :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` denotes the + matrix elements of the operator :math:`\hat{v}` + + .. math:: + + \langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle = + \int dr_1 \int dr_2 ~ \phi_\alpha^*(r_1) \phi_\beta^*(r_2) ~\hat{v}(r_1, r_2)~ + \phi_\gamma(r_2) \phi_\delta(r_1). + + If an active space is defined (see :func:`~.active_space`), the summation indices + run over the active orbitals and the contribution due to core orbitals is computed as + + .. math:: + + && \hat{V}_\mathrm{core} = v_\mathrm{core} + + \sum_{\alpha, \beta \in \mathrm{active}} \sum_{i \in \mathrm{core}} + (2 \langle i, \alpha \vert \hat{v} \vert \beta, i \rangle - + \langle i, \alpha \vert \hat{v} \vert i, \beta \rangle)~ + [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} + + \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}] \\ + && v_\mathrm{core} = \sum_{\alpha,\beta \in \mathrm{core}} + [2 \langle \alpha, \beta \vert \hat{v} \vert \beta, \alpha \rangle - + \langle \alpha, \beta \vert \hat{v} \vert \alpha, \beta \rangle]. + + Args: + matrix_elements (array[float]): 4D NumPy array with the matrix elements + :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` + core (list): indices of core orbitals, i.e., the orbitals that are + not correlated in the many-body wave function + active (list): indices of active orbitals, i.e., the orbitals used to + build the correlated many-body wave function + cutoff (float): Cutoff value for including matrix elements. The + matrix elements with absolute value less than ``cutoff`` are neglected. + + Returns: + FermionOperator: an instance of OpenFermion's ``FermionOperator`` representing the + two-particle operator :math:`\hat{V}`. + + **Example** + + >>> import numpy as np + >>> matrix_elements = np.array([[[[ 6.82389533e-01, -1.45716772e-16], + ... [-2.77555756e-17, 1.79000576e-01]], + ... [[-2.77555756e-17, 1.79000576e-16], + ... [ 6.70732778e-01, 0.00000000e+00]]], + ... [[[-1.45716772e-16, 6.70732778e-16], + ... [ 1.79000576e-01, -8.32667268e-17]], + ... [[ 1.79000576e-16, -8.32667268e-17], + ... [ 0.00000000e+00, 7.05105632e-01]]]]) + >>> v_op = two_particle(matrix_elements) + >>> print(v_op) + 0.3411947665 [0^ 0^ 0 0] + + 0.089500288 [0^ 0^ 2 2] + + 0.3411947665 [0^ 1^ 1 0] + + 0.089500288 [0^ 1^ 3 2] + + 0.335366389 [0^ 2^ 2 0] + + 0.335366389 [0^ 3^ 3 0] + + 0.3411947665 [1^ 0^ 0 1] + + 0.089500288 [1^ 0^ 2 3] + + 0.3411947665 [1^ 1^ 1 1] + + 0.089500288 [1^ 1^ 3 3] + + 0.335366389 [1^ 2^ 2 1] + + 0.335366389 [1^ 3^ 3 1] + + 0.089500288 [2^ 0^ 2 0] + + 0.089500288 [2^ 1^ 3 0] + + 0.352552816 [2^ 2^ 2 2] + + 0.352552816 [2^ 3^ 3 2] + + 0.089500288 [3^ 0^ 2 1] + + 0.089500288 [3^ 1^ 3 1] + + 0.352552816 [3^ 2^ 2 3] + + 0.352552816 [3^ 3^ 3 3] + """ + openfermion, _ = _import_of() + + orbitals = matrix_elements.shape[0] + + if matrix_elements.ndim != 4: + raise ValueError( + f"'matrix_elements' must be a 4D array; got 'matrix_elements.ndim = ' {matrix_elements.ndim}" + ) + + if not core: + v_core = 0 + else: + if any(i > orbitals - 1 or i < 0 for i in core): + raise ValueError( + f"Indices of core orbitals must be between 0 and {orbitals - 1}; got core = {core}" + ) + + # Compute the contribution of core orbitals + v_core = sum( + [ + 2 * matrix_elements[alpha, beta, beta, alpha] + - matrix_elements[alpha, beta, alpha, beta] + for alpha in core + for beta in core + ] + ) + + if active is None: + if not core: + active = list(range(orbitals)) + else: + active = [i for i in range(orbitals) if i not in core] + + if any(i > orbitals - 1 or i < 0 for i in active): + raise ValueError( + f"Indices of active orbitals must be between 0 and {orbitals - 1}; got active = {active}" + ) + + # Indices of the matrix elements with absolute values >= cutoff + indices = np.nonzero(np.abs(matrix_elements) >= cutoff) + + # Single out the indices of active orbitals + num_indices = len(indices[0]) + quads = [ + [indices[0][i], indices[1][i], indices[2][i], indices[3][i]] + for i in range(num_indices) + if all(indices[j][i] in active for j in range(len(indices))) + ] + + # Build the FermionOperator representing V + v_op = openfermion.ops.FermionOperator("") * v_core + + # add renormalized (due to core orbitals) "one-particle" operators + if core: + for alpha in active: + for beta in active: + element = 2 * np.sum( + matrix_elements[np.array(core), alpha, beta, np.array(core)] + ) - np.sum(matrix_elements[np.array(core), alpha, np.array(core), beta]) + + # up-up term + a = 2 * active.index(alpha) + b = 2 * active.index(beta) + v_op += openfermion.ops.FermionOperator(((a, 1), (b, 0)), element) + + # down-down term + v_op += openfermion.ops.FermionOperator(((a + 1, 1), (b + 1, 0)), element) + + # add two-particle operators + for quad in quads: + alpha, beta, gamma, delta = quad + element = matrix_elements[alpha, beta, gamma, delta] + + # up-up-up-up term + a = 2 * active.index(alpha) + b = 2 * active.index(beta) + g = 2 * active.index(gamma) + d = 2 * active.index(delta) + v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) + + # up-down-down-up term + b = 2 * active.index(beta) + 1 + g = 2 * active.index(gamma) + 1 + v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) + + # down-up-up-down term + a = 2 * active.index(alpha) + 1 + b = 2 * active.index(beta) + g = 2 * active.index(gamma) + d = 2 * active.index(delta) + 1 + v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) + + # down-down-down-down term + b = 2 * active.index(beta) + 1 + g = 2 * active.index(gamma) + 1 + v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) + + return v_op + + +def dipole_of( + symbols, + coordinates, + name="molecule", + charge=0, + mult=1, + basis="sto-3g", + package="pyscf", + core=None, + active=None, + mapping="jordan_wigner", + cutoff=1.0e-12, + outpath=".", + wires=None, +): + r"""Computes the electric dipole moment operator in the Pauli basis. + + The second quantized dipole moment operator :math:`\hat{D}` of a molecule is given by + + .. math:: + + \hat{D} = -\sum_{\alpha, \beta} \langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle + [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} + + \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}] + \hat{D}_\mathrm{n}. + + In the equation above, the indices :math:`\alpha, \beta` run over the basis of Hartree-Fock + molecular orbitals and the operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the + electron creation and annihilation operators, respectively. The matrix elements of the + position operator :math:`\hat{{\bf r}}` are computed as + + .. math:: + + \langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle = \sum_{i, j} + C_{\alpha i}^*C_{\beta j} \langle i \vert \hat{{\bf r}} \vert j \rangle, + + where :math:`\vert i \rangle` is the wave function of the atomic orbital, + :math:`C_{\alpha i}` are the coefficients defining the molecular orbitals, + and :math:`\langle i \vert \hat{{\bf r}} \vert j \rangle` + is the representation of operator :math:`\hat{{\bf r}}` in the atomic basis. + + The contribution of the nuclei to the dipole operator is given by + + .. math:: + + \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i \hat{I}, + + + where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the + nuclear coordinates of the :math:`i`-th atom of the molecule. + + Args: + symbols (list[str]): symbols of the atomic species in the molecule + coordinates (array[float]): 1D array with the atomic positions in Cartesian + coordinates. The coordinates must be given in atomic units and the size of the array + should be ``3*N`` where ``N`` is the number of atoms. + name (str): name of the molecule + charge (int): charge of the molecule + mult (int): spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1` of the + Hartree-Fock (HF) state based on the number of unpaired electrons occupying the + HF orbitals + basis (str): Atomic basis set used to represent the molecular orbitals. Basis set + availability per element can be found + `here `_ + package (str): quantum chemistry package (pyscf) used to solve the + mean field electronic structure problem + core (list): indices of core orbitals + active (list): indices of active orbitals + mapping (str): transformation (``'jordan_wigner'`` or ``'bravyi_kitaev'``) used to + map the fermionic operator to the Pauli basis + cutoff (float): Cutoff value for including the matrix elements + :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements + with absolute value less than ``cutoff`` are neglected. + outpath (str): path to the directory containing output files + wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator + to an observable measurable in a PennyLane ansatz. + For types Wires/list/tuple, each item in the iterable represents a wire label + corresponding to the qubit number equal to its index. + For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted. + If None, will use identity map (e.g. 0->0, 1->1, ...). + + Returns: + list[pennylane.Hamiltonian]: the qubit observables corresponding to the components + :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator in + atomic units. + + **Example** + + >>> symbols = ["H", "H", "H"] + >>> coordinates = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0, 1.855, 0.002, 0.0]) + >>> dipole_obs = dipole_of(symbols, coordinates, charge=1) + >>> print([(h.wires) for h in dipole_obs]) + [, , ] + + >>> dipole_obs[0] # x-component of D + ( + 0.4781123173263876 * Z(0) + + 0.4781123173263876 * Z(1) + + -0.3913638489489803 * (Y(0) @ Z(1) @ Y(2)) + + -0.3913638489489803 * (X(0) @ Z(1) @ X(2)) + + -0.3913638489489803 * (Y(1) @ Z(2) @ Y(3)) + + -0.3913638489489803 * (X(1) @ Z(2) @ X(3)) + + 0.2661114704527088 * (Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4)) + + 0.2661114704527088 * (X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4)) + + 0.2661114704527088 * (Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5)) + + 0.2661114704527088 * (X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5)) + + 0.7144779061810713 * Z(2) + + 0.7144779061810713 * Z(3) + + -0.11734958781031017 * (Y(2) @ Z(3) @ Y(4)) + + -0.11734958781031017 * (X(2) @ Z(3) @ X(4)) + + -0.11734958781031017 * (Y(3) @ Z(4) @ Y(5)) + + -0.11734958781031017 * (X(3) @ Z(4) @ X(5)) + + 0.24190977644645698 * Z(4) + + 0.24190977644645698 * Z(5) + ) + """ + openfermion, _ = _import_of() + + if mult != 1: + raise ValueError( + f"Currently, this functionality is constrained to Hartree-Fock states with spin multiplicity = 1;" + f" got multiplicity 2S+1 = {mult}" + ) + + for i in symbols: + if i not in atomic_numbers: + raise ValueError( + f"Currently, only first- or second-row elements of the periodic table are supported;" + f" got element {i}" + ) + + hf_file = qml.qchem.meanfield(symbols, coordinates, name, charge, mult, basis, package, outpath) + + hf = openfermion.MolecularData(filename=hf_file.strip()) + + # Load dipole matrix elements in the atomic basis + # pylint: disable=import-outside-toplevel + from pyscf import gto + + mol = gto.M( + atom=hf.geometry, basis=hf.basis, charge=hf.charge, spin=0.5 * (hf.multiplicity - 1) + ) + dip_ao = mol.intor_symmetric("int1e_r", comp=3).real + + # Transform dipole matrix elements to the MO basis + n_orbs = hf.n_orbitals + c_hf = hf.canonical_orbitals + + dip_mo = np.zeros((3, n_orbs, n_orbs)) + for comp in range(3): + for alpha in range(n_orbs): + for beta in range(alpha + 1): + dip_mo[comp, alpha, beta] = c_hf[:, alpha] @ dip_ao[comp] @ c_hf[:, beta] + + dip_mo[comp] += dip_mo[comp].T - np.diag(np.diag(dip_mo[comp])) + + # Compute the nuclear contribution + dip_n = np.zeros(3) + for comp in range(3): + for i, symb in enumerate(symbols): + dip_n[comp] += atomic_numbers[symb] * coordinates[3 * i + comp] + + # Build the observable + dip = [] + for i in range(3): + fermion_obs = one_particle(dip_mo[i], core=core, active=active, cutoff=cutoff) + dip.append(observable([-fermion_obs], init_term=dip_n[i], mapping=mapping, wires=wires)) + + return dip + + +def molecular_dipole( + molecule, + method="dhf", + active_electrons=None, + active_orbitals=None, + mapping="jordan_wigner", + outpath=".", + wires=None, + args=None, + cutoff=1.0e-16, +): # pylint:disable=too-many-arguments, too-many-statements + r"""Generate the dipole moment operator for a molecule in the Pauli basis. + + The dipole operator in the second-quantized form is + + .. math:: + + \hat{D} = -\sum_{pq} d_{pq} [\hat{c}_{p\uparrow}^\dagger \hat{c}_{q\uparrow} + + \hat{c}_{p\downarrow}^\dagger \hat{c}_{q\downarrow}] - + \hat{D}_\mathrm{c} + \hat{D}_\mathrm{n}, + + where the matrix elements :math:`d_{pq}` are given by the integral of the position operator + :math:`\hat{{\bf r}}` over molecular orbitals :math:`\phi` + + .. math:: + + d_{pq} = \int \phi_p^*(r) \hat{{\bf r}} \phi_q(r) dr, + + and :math:`\hat{c}^{\dagger}` and :math:`\hat{c}` are the creation and annihilation operators, + respectively. The contribution of the core orbitals and nuclei are denoted by + :math:`\hat{D}_\mathrm{c}` and :math:`\hat{D}_\mathrm{n}`, respectively, which are computed as + + .. math:: + \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii} \quad \text{and} \quad + \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i, + + where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the + nuclear coordinates of the :math:`i`-th atom of the molecule. + + The fermionic dipole operator is then transformed to the qubit basis, which gives + + .. math:: + + \hat{D} = \sum_{j} c_j P_j, + + where :math:`c_j` is a numerical coefficient and :math:`P_j` is a tensor product of + single-qubit Pauli operators :math:`X, Y, Z, I`. The qubit observables corresponding + to the components :math:`\hat{D}_x`, :math:`\hat{D}_y`, and :math:`\hat{D}_z` of the + dipole operator are then computed separately. + + Args: + 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, 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 + are considered to be active. + mapping (str): Transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian. + Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. + outpath (str): Path to the directory containing output files + wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator to + an observable measurable in a Pennylane ansatz. + For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label + 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. + args (array[array[float]]): Initial values of the differentiable parameters + cutoff (float): Cutoff value for including the matrix elements + :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements + with absolute value less than ``cutoff`` are neglected. + + Returns: + list[pennylane.Hamiltonian]: The qubit observables corresponding to the components + :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator. + + **Example** + + >>> symbols = ["H", "H", "H"] + >>> coordinates = np.array([[0.028, 0.054, 0.0], [0.986, 1.610, 0.0], [1.855, 0.002, 0.0]]) + >>> mol = qml.qchem.Molecule(symbols, coordinates, charge=1) + >>> dipole_obs = qml.qchem.molecular_dipole(mol, method="openfermion") + >>> dipole_obs[0] # x-component of D + ( + 0.4781123173263876 * Z(0) + + 0.4781123173263876 * Z(1) + + -0.3913638489489803 * (Y(0) @ Z(1) @ Y(2)) + + -0.3913638489489803 * (X(0) @ Z(1) @ X(2)) + + -0.3913638489489803 * (Y(1) @ Z(2) @ Y(3)) + + -0.3913638489489803 * (X(1) @ Z(2) @ X(3)) + + 0.2661114704527088 * (Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4)) + + 0.2661114704527088 * (X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4)) + + 0.2661114704527088 * (Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5)) + + 0.2661114704527088 * (X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5)) + + 0.7144779061810713 * Z(2) + + 0.7144779061810713 * Z(3) + + -0.11734958781031017 * (Y(2) @ Z(3) @ Y(4)) + + -0.11734958781031017 * (X(2) @ Z(3) @ X(4)) + + -0.11734958781031017 * (Y(3) @ Z(4) @ Y(5)) + + -0.11734958781031017 * (X(3) @ Z(4) @ X(5)) + + 0.24190977644645698 * Z(4) + + 0.24190977644645698 * Z(5) + ) + """ + + if method not in ["dhf", "openfermion"]: + raise ValueError("Only 'dhf', 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'." + ) + + symbols = molecule.symbols + coordinates = molecule.coordinates + + if qml.math.shape(coordinates)[0] == len(symbols) * 3: + geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3)) + geometry_hf = coordinates + elif len(coordinates) == len(symbols): + geometry_dhf = qml.numpy.array(coordinates) + geometry_hf = coordinates.flatten() + + if molecule.mult != 1: + raise ValueError( + "Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule." + ) + + core, active = qml.qchem.active_space( + molecule.n_electrons, molecule.n_orbitals, molecule.mult, active_electrons, active_orbitals + ) + + if method == "dhf": + + if args is None and isinstance(geometry_dhf, qml.numpy.tensor): + geometry_dhf.requires_grad = False + mol = qml.qchem.Molecule( + symbols, + geometry_dhf, + charge=molecule.charge, + mult=molecule.mult, + basis_name=molecule.basis_name, + load_data=molecule.load_data, + alpha=molecule.alpha, + coeff=molecule.coeff, + ) + + requires_grad = args is not None + dip = ( + qml.qchem.dipole_moment(mol, cutoff=cutoff, core=core, active=active, mapping=mapping)( + *(args or []) + ) + if requires_grad + else qml.qchem.dipole_moment( + mol, cutoff=cutoff, core=core, active=active, mapping=mapping + )() + ) + if wires: + wires_new = qml.qchem.convert._process_wires(wires) + wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels))) + dip = [qml.map_wires(op, wires_map) for op in dip] + + return dip + + dip = qml.qchem.dipole_of( + symbols, + geometry_hf, + molecule.name, + molecule.charge, + molecule.mult, + molecule.basis_name, + package="pyscf", + core=core, + active=active, + mapping=mapping, + cutoff=cutoff, + outpath=outpath, + wires=wires, + ) + + return dip + + +def meanfield( + symbols, + coordinates, + name="molecule", + charge=0, + mult=1, + basis="sto-3g", + package="pyscf", + outpath=".", +): # pylint: disable=too-many-arguments + r"""Generates a file from which the mean field electronic structure + of the molecule can be retrieved. + + This function uses OpenFermion-PySCF plugins to + perform the Hartree-Fock (HF) calculation for the polyatomic system using the quantum + chemistry packages ``PySCF``. The mean field electronic + structure is saved in an hdf5-formatted file. + + The charge of the molecule can be given to simulate cationic/anionic systems. + Also, the spin multiplicity can be input to determine the number of unpaired electrons + occupying the HF orbitals as illustrated in the figure below. + + | + + .. figure:: ../../_static/qchem/hf_references.png + :align: center + :width: 50% + + | + + Args: + symbols (list[str]): symbols of the atomic species in the molecule + coordinates (array[float]): 1D array with the atomic positions in Cartesian + coordinates. The coordinates must be given in atomic units and the size of the array + should be ``3*N`` where ``N`` is the number of atoms. + name (str): molecule label + charge (int): net charge of the system + mult (int): Spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1` for + :math:`N_\mathrm{unpaired}` unpaired electrons occupying the HF orbitals. + Possible values for ``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 HF orbitals. Basis set + availability per element can be found + `here `_ + package (str): Quantum chemistry package used to solve the Hartree-Fock equations. + outpath (str): path to output directory + + Returns: + str: absolute path to the file containing the mean field electronic structure + + **Example** + + >>> symbols, coordinates = (['H', 'H'], np.array([0., 0., -0.66140414, 0., 0., 0.66140414])) + >>> meanfield(symbols, coordinates, name="h2") + ./h2_pyscf_sto-3g + """ + openfermion, openfermionpyscf = _import_of() + + if coordinates.size != 3 * len(symbols): + raise ValueError( + f"The size of the array 'coordinates' has to be 3*len(symbols) = {3 * len(symbols)};" + f" got 'coordinates.size' = {coordinates.size}" + ) + + package = package.strip().lower() + + if package not in "pyscf": + error_message = ( + f"Integration with quantum chemistry package '{package}' is not available. \n Please set" + f" 'package' to 'pyscf'." + ) + raise TypeError(error_message) + + filename = name + "_" + package.lower() + "_" + basis.strip() + path_to_file = os.path.join(outpath.strip(), filename) + + geometry = [ + [symbol, tuple(np.array(coordinates)[3 * i : 3 * i + 3] * bohr_angs)] + for i, symbol in enumerate(symbols) + ] + + molecule = openfermion.MolecularData(geometry, basis, mult, charge, filename=path_to_file) + + if package == "pyscf": + # pylint: disable=import-outside-toplevel + from openfermionpyscf import run_pyscf + + run_pyscf(molecule, run_scf=1, verbose=0) + + return path_to_file + + +def decompose(hf_file, mapping="jordan_wigner", core=None, active=None): + r"""Decomposes the molecular Hamiltonian into a linear combination of Pauli operators using + OpenFermion tools. + + This function uses OpenFermion functions to build the second-quantized electronic Hamiltonian + of the molecule and map it to the Pauli basis using the Jordan-Wigner, Parity or Bravyi-Kitaev + transformation. + + Args: + hf_file (str): absolute path to the hdf5-formatted file with the + Hartree-Fock electronic structure + mapping (str): Specifies the transformation to map the fermionic Hamiltonian to the + Pauli basis. Input values can be ``'jordan_wigner'``, ``'parity'``, or ``'bravyi_kitaev'``. + core (list): indices of core orbitals, i.e., the orbitals that are + not correlated in the many-body wave function + active (list): indices of active orbitals, i.e., the orbitals used to + build the correlated many-body wave function + + Returns: + QubitOperator: an instance of OpenFermion's ``QubitOperator`` + + **Example** + + >>> decompose('./pyscf/sto-3g/h2', mapping='bravyi_kitaev') + (-0.04207897696293986+0j) [] + (0.04475014401986122+0j) [X0 Z1 X2] + + (0.04475014401986122+0j) [X0 Z1 X2 Z3] +(0.04475014401986122+0j) [Y0 Z1 Y2] + + (0.04475014401986122+0j) [Y0 Z1 Y2 Z3] +(0.17771287459806262+0j) [Z0] + + (0.17771287459806265+0j) [Z0 Z1] +(0.1676831945625423+0j) [Z0 Z1 Z2] + + (0.1676831945625423+0j) [Z0 Z1 Z2 Z3] +(0.12293305054268105+0j) [Z0 Z2] + + (0.12293305054268105+0j) [Z0 Z2 Z3] +(0.1705973832722409+0j) [Z1] + + (-0.2427428049645989+0j) [Z1 Z2 Z3] +(0.1762764080276107+0j) [Z1 Z3] + + (-0.2427428049645989+0j) [Z2] + """ + openfermion, _ = _import_of() + + # loading HF data from the hdf5 file + molecule = openfermion.MolecularData(filename=hf_file.strip()) + + # getting the terms entering the second-quantized Hamiltonian + terms_molecular_hamiltonian = molecule.get_molecular_hamiltonian( + occupied_indices=core, active_indices=active + ) + + # generating the fermionic Hamiltonian + fermionic_hamiltonian = openfermion.transforms.get_fermion_operator(terms_molecular_hamiltonian) + + mapping = mapping.strip().lower() + + # 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 _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 + + +def _openfermion_hamiltonian( + symbols, + coordinates, + name, + charge=0, + mult=1, + basis="sto-3g", + active_electrons=None, + active_orbitals=None, + mapping="jordan_wigner", + outpath=".", + wires=None, + convert_tol=1e012, +): + r"""Compute openfermion Hamiltonian.""" + openfermion, _ = _import_of() + + hf_file = meanfield(symbols, coordinates, name, charge, mult, basis, "pyscf", outpath) + + molecule = openfermion.MolecularData(filename=hf_file) + + core, active = qml.qchem.active_space( + molecule.n_electrons, molecule.n_orbitals, mult, active_electrons, active_orbitals + ) + + h_of, qubits = (decompose(hf_file, mapping, core, active), 2 * len(active)) + + h_pl = qml.qchem.convert.import_operator(h_of, wires=wires, tol=convert_tol) + + return h_pl From 607be7a84a472e444beebb7de0a4b6f683f6c3ee Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Thu, 6 Jun 2024 12:45:49 -0400 Subject: [PATCH 26/28] Resolve merge conflicts --- pennylane/qchem/openfermion_obs.py | 1434 ---------------------------- 1 file changed, 1434 deletions(-) delete mode 100644 pennylane/qchem/openfermion_obs.py diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py deleted file mode 100644 index 70a2c078d63..00000000000 --- a/pennylane/qchem/openfermion_obs.py +++ /dev/null @@ -1,1434 +0,0 @@ -# Copyright 2018-2020 Xanadu Quantum Technologies Inc. - -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at - -# http://www.apache.org/licenses/LICENSE-2.0 - -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -"""This module contains functions to construct many-body observables with ``OpenFermion-PySCF``. -""" -# 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 - - -def _import_of(): - """Import openfermion and openfermionpyscf.""" - try: - # pylint: disable=import-outside-toplevel, unused-import, multiple-imports - import openfermion - import openfermionpyscf - except ImportError as Error: - raise ImportError( - "This feature requires openfermionpyscf. " - "It can be installed with: pip install openfermionpyscf." - ) from Error - - 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. - - The second-quantized operator of the fermionic many-body system can combine one-particle - and two-particle operators as in the case of electronic Hamiltonians :math:`\hat{H}`: - - .. math:: - - \hat{H} = \sum_{\alpha, \beta} \langle \alpha \vert \hat{t}^{(1)} + - \cdots + \hat{t}^{(n)} \vert \beta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta - + \frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta} - \langle \alpha, \beta \vert \hat{v}^{(1)} + \cdots + \hat{v}^{(n)} - \vert \gamma, \delta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta^\dagger - \hat{c}_\gamma \hat{c}_\delta - - In the latter equations the indices :math:`\alpha, \beta, \gamma, \delta` run over the - basis of single-particle states. The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` - are the particle creation and annihilation operators, respectively. - :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` denotes the matrix element of - the single-particle operator :math:`\hat{t}` entering the observable. For example, - in electronic structure calculations, this is the case for the kinetic energy operator, - the nuclei Coulomb potential, or any other external fields included in the Hamiltonian. - On the other hand, :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` - denotes the matrix element of the two-particle operator :math:`\hat{v}`, for example, the - Coulomb interaction between the electrons. - - - The observable is built by adding the operators - :math:`\sum_{\alpha, \beta} t_{\alpha\beta}^{(i)} - \hat{c}_\alpha^\dagger \hat{c}_\beta` and - :math:`\frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta} - v_{\alpha\beta\gamma\delta}^{(i)} - \hat{c}_\alpha^\dagger \hat{c}_\beta^\dagger \hat{c}_\gamma \hat{c}_\delta`. - - - Second-quantized operators contributing to the - many-body observable must be represented using the `FermionOperator - `_ data structure as implemented in OpenFermion. - See the functions :func:`~.one_particle` and :func:`~.two_particle` to build the - FermionOperator representations of one-particle and two-particle operators. - - - The function uses tools of `OpenFermion `_ - to map the resulting fermionic Hamiltonian to the basis of Pauli matrices via the - Jordan-Wigner or Bravyi-Kitaev transformation. Finally, the qubit operator is converted - to a PennyLane observable by the function :func:`~.convert_observable`. - - Args: - fermion_ops (list[FermionOperator]): list containing the FermionOperator data structures - representing the one-particle and/or two-particle operators entering the many-body - observable - init_term (float): Any quantity required to initialize the many-body observable. For - example, this can be used to pass the nuclear-nuclear repulsion energy :math:`V_{nn}` - which is typically included in the electronic Hamiltonian of molecules. - mapping (str): Specifies the fermion-to-qubit mapping. Input values can - be ``'jordan_wigner'``, ``'parity'``, or ``'bravyi_kitaev'``. - wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator - to an observable measurable in a PennyLane ansatz. - For types Wires/list/tuple, each item in the iterable represents a wire label - corresponding to the qubit number equal to its index. - For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted. - If None, will use identity map (e.g. 0->0, 1->1, ...). - - Returns: - pennylane.Hamiltonian: the fermionic-to-qubit transformed observable - - **Example** - - >>> t = FermionOperator("0^ 0", 0.5) + FermionOperator("1^ 1", 0.25) - >>> v = FermionOperator("1^ 0^ 0 1", -0.15) + FermionOperator("2^ 0^ 2 0", 0.3) - >>> observable([t, v], mapping="jordan_wigner") - ( - 0.2625 * I(0) - + -0.1375 * Z(0) - + -0.0875 * Z(1) - + -0.0375 * (Z(0) @ Z(1)) - + 0.075 * Z(2) - + -0.075 * (Z(0) @ Z(2)) - ) - """ - openfermion, _ = _import_of() - - if mapping.strip().lower() not in ("jordan_wigner", "parity", "bravyi_kitaev"): - raise TypeError( - f"The '{mapping}' transformation is not available. \n " - f"Please set 'mapping' to 'jordan_wigner', 'parity', or 'bravyi_kitaev'." - ) - - # Initialize the FermionOperator - mb_obs = openfermion.ops.FermionOperator("") * init_term - for ops in fermion_ops: - if not isinstance(ops, openfermion.ops.FermionOperator): - raise TypeError( - f"Elements in the lists are expected to be of type 'FermionOperator'; got {type(ops)}" - ) - mb_obs += ops - - # Map the fermionic operator to a qubit operator - if mapping.strip().lower() == "bravyi_kitaev": - return qml.qchem.convert.import_operator( - openfermion.transforms.bravyi_kitaev(mb_obs), wires=wires - ) - - if mapping == "parity": - qubits = openfermion.count_qubits(mb_obs) - if qubits == 0: - return 0.0 * qml.I(0) - binary_code = openfermion.parity_code(qubits) - return qml.qchem.convert.import_operator( - openfermion.transforms.binary_code_transform(mb_obs, binary_code), wires=wires - ) - - return qml.qchem.convert.import_operator( - openfermion.transforms.jordan_wigner(mb_obs), wires=wires - ) - - -def one_particle(matrix_elements, core=None, active=None, cutoff=1.0e-12): - r"""Generates the `FermionOperator `_ representing a given one-particle operator - required to build many-body qubit observables. - - Second quantized one-particle operators are expanded in the basis of single-particle - states as - - .. math:: - - \hat{T} = \sum_{\alpha, \beta} \langle \alpha \vert \hat{t} \vert \beta \rangle - [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} + - \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}]. - - In the equation above the indices :math:`\alpha, \beta` run over the basis of spatial - orbitals :math:`\phi_\alpha(r)`. Since the operator :math:`\hat{t}` acts only on the - spatial coordinates, the spin quantum numbers are indicated explicitly with the up/down arrows. - The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation - and annihilation operators, respectively, and - :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` denotes the matrix elements of - the operator :math:`\hat{t}` - - .. math:: - - \langle \alpha \vert \hat{t} \vert \beta \rangle = \int dr ~ \phi_\alpha^*(r) - \hat{t}(r) \phi_\beta(r). - - If an active space is defined (see :func:`~.active_space`), the summation indices - run over the active orbitals and the contribution due to core orbitals is computed as - :math:`t_\mathrm{core} = 2 \sum_{\alpha\in \mathrm{core}} - \langle \alpha \vert \hat{t} \vert \alpha \rangle`. - - Args: - matrix_elements (array[float]): 2D NumPy array with the matrix elements - :math:`\langle \alpha \vert \hat{t} \vert \beta \rangle` - core (list): indices of core orbitals, i.e., the orbitals that are - not correlated in the many-body wave function - active (list): indices of active orbitals, i.e., the orbitals used to - build the correlated many-body wave function - cutoff (float): Cutoff value for including matrix elements. The - matrix elements with absolute value less than ``cutoff`` are neglected. - - Returns: - FermionOperator: an instance of OpenFermion's ``FermionOperator`` representing the - one-particle operator :math:`\hat{T}`. - - **Example** - - >>> import numpy as np - >>> matrix_elements = np.array([[-1.27785301e+00, 0.00000000e+00], - ... [ 1.52655666e-16, -4.48299696e-01]]) - >>> t_op = one_particle(matrix_elements) - >>> print(t_op) - -1.277853006156875 [0^ 0] + - -1.277853006156875 [1^ 1] + - -0.44829969610163756 [2^ 2] + - -0.44829969610163756 [3^ 3] - """ - openfermion, _ = _import_of() - - orbitals = matrix_elements.shape[0] - - if matrix_elements.ndim != 2: - raise ValueError( - f"'matrix_elements' must be a 2D array; got matrix_elements.ndim = {matrix_elements.ndim}" - ) - - if not core: - t_core = 0 - else: - if any(i > orbitals - 1 or i < 0 for i in core): - raise ValueError( - f"Indices of core orbitals must be between 0 and {orbitals - 1}; got core = {core}" - ) - - # Compute contribution due to core orbitals - t_core = 2 * np.sum(matrix_elements[np.array(core), np.array(core)]) - - if active is None: - if not core: - active = list(range(orbitals)) - else: - active = [i for i in range(orbitals) if i not in core] - - if any(i > orbitals - 1 or i < 0 for i in active): - raise ValueError( - f"Indices of active orbitals must be between 0 and {orbitals - 1}; got active = {active}" - ) - - # Indices of the matrix elements with absolute values >= cutoff - indices = np.nonzero(np.abs(matrix_elements) >= cutoff) - - # Single out the indices of active orbitals - num_indices = len(indices[0]) - pairs = [ - [indices[0][i], indices[1][i]] - for i in range(num_indices) - if all(indices[j][i] in active for j in range(len(indices))) - ] - - # Build the FermionOperator representing T - t_op = openfermion.ops.FermionOperator("") * t_core - for pair in pairs: - alpha, beta = pair - element = matrix_elements[alpha, beta] - - # spin-up term - a = 2 * active.index(alpha) - b = 2 * active.index(beta) - t_op += openfermion.ops.FermionOperator(((a, 1), (b, 0)), element) - - # spin-down term - t_op += openfermion.ops.FermionOperator(((a + 1, 1), (b + 1, 0)), element) - - return t_op - - -def two_particle(matrix_elements, core=None, active=None, cutoff=1.0e-12): - r"""Generates the `FermionOperator `_ representing a given two-particle operator - required to build many-body qubit observables. - - Second quantized two-particle operators are expanded in the basis of single-particle - states as - - .. math:: - - \hat{V} = \frac{1}{2} \sum_{\alpha, \beta, \gamma, \delta} - \langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle - ~ &[& \hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow}^\dagger - \hat{c}_{\gamma\uparrow} \hat{c}_{\delta\uparrow} + \hat{c}_{\alpha\uparrow}^\dagger - \hat{c}_{\beta\downarrow}^\dagger \hat{c}_{\gamma\downarrow} \hat{c}_{\delta\uparrow} \\ - &+& \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\uparrow}^\dagger - \hat{c}_{\gamma\uparrow} \hat{c}_{\delta\downarrow} + \hat{c}_{\alpha\downarrow}^\dagger - \hat{c}_{\beta\downarrow}^\dagger \hat{c}_{\gamma\downarrow} \hat{c}_{\delta\downarrow}~]. - - In the equation above the indices :math:`\alpha, \beta, \gamma, \delta` run over the basis - of spatial orbitals :math:`\phi_\alpha(r)`. Since the operator :math:`v` acts only on the - spatial coordinates the spin quantum numbers are indicated explicitly with the up/down arrows. - The operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the particle creation and - annihilation operators, respectively, and - :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` denotes the - matrix elements of the operator :math:`\hat{v}` - - .. math:: - - \langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle = - \int dr_1 \int dr_2 ~ \phi_\alpha^*(r_1) \phi_\beta^*(r_2) ~\hat{v}(r_1, r_2)~ - \phi_\gamma(r_2) \phi_\delta(r_1). - - If an active space is defined (see :func:`~.active_space`), the summation indices - run over the active orbitals and the contribution due to core orbitals is computed as - - .. math:: - - && \hat{V}_\mathrm{core} = v_\mathrm{core} + - \sum_{\alpha, \beta \in \mathrm{active}} \sum_{i \in \mathrm{core}} - (2 \langle i, \alpha \vert \hat{v} \vert \beta, i \rangle - - \langle i, \alpha \vert \hat{v} \vert i, \beta \rangle)~ - [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} + - \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}] \\ - && v_\mathrm{core} = \sum_{\alpha,\beta \in \mathrm{core}} - [2 \langle \alpha, \beta \vert \hat{v} \vert \beta, \alpha \rangle - - \langle \alpha, \beta \vert \hat{v} \vert \alpha, \beta \rangle]. - - Args: - matrix_elements (array[float]): 4D NumPy array with the matrix elements - :math:`\langle \alpha, \beta \vert \hat{v} \vert \gamma, \delta \rangle` - core (list): indices of core orbitals, i.e., the orbitals that are - not correlated in the many-body wave function - active (list): indices of active orbitals, i.e., the orbitals used to - build the correlated many-body wave function - cutoff (float): Cutoff value for including matrix elements. The - matrix elements with absolute value less than ``cutoff`` are neglected. - - Returns: - FermionOperator: an instance of OpenFermion's ``FermionOperator`` representing the - two-particle operator :math:`\hat{V}`. - - **Example** - - >>> import numpy as np - >>> matrix_elements = np.array([[[[ 6.82389533e-01, -1.45716772e-16], - ... [-2.77555756e-17, 1.79000576e-01]], - ... [[-2.77555756e-17, 1.79000576e-16], - ... [ 6.70732778e-01, 0.00000000e+00]]], - ... [[[-1.45716772e-16, 6.70732778e-16], - ... [ 1.79000576e-01, -8.32667268e-17]], - ... [[ 1.79000576e-16, -8.32667268e-17], - ... [ 0.00000000e+00, 7.05105632e-01]]]]) - >>> v_op = two_particle(matrix_elements) - >>> print(v_op) - 0.3411947665 [0^ 0^ 0 0] + - 0.089500288 [0^ 0^ 2 2] + - 0.3411947665 [0^ 1^ 1 0] + - 0.089500288 [0^ 1^ 3 2] + - 0.335366389 [0^ 2^ 2 0] + - 0.335366389 [0^ 3^ 3 0] + - 0.3411947665 [1^ 0^ 0 1] + - 0.089500288 [1^ 0^ 2 3] + - 0.3411947665 [1^ 1^ 1 1] + - 0.089500288 [1^ 1^ 3 3] + - 0.335366389 [1^ 2^ 2 1] + - 0.335366389 [1^ 3^ 3 1] + - 0.089500288 [2^ 0^ 2 0] + - 0.089500288 [2^ 1^ 3 0] + - 0.352552816 [2^ 2^ 2 2] + - 0.352552816 [2^ 3^ 3 2] + - 0.089500288 [3^ 0^ 2 1] + - 0.089500288 [3^ 1^ 3 1] + - 0.352552816 [3^ 2^ 2 3] + - 0.352552816 [3^ 3^ 3 3] - """ - openfermion, _ = _import_of() - - orbitals = matrix_elements.shape[0] - - if matrix_elements.ndim != 4: - raise ValueError( - f"'matrix_elements' must be a 4D array; got 'matrix_elements.ndim = ' {matrix_elements.ndim}" - ) - - if not core: - v_core = 0 - else: - if any(i > orbitals - 1 or i < 0 for i in core): - raise ValueError( - f"Indices of core orbitals must be between 0 and {orbitals - 1}; got core = {core}" - ) - - # Compute the contribution of core orbitals - v_core = sum( - [ - 2 * matrix_elements[alpha, beta, beta, alpha] - - matrix_elements[alpha, beta, alpha, beta] - for alpha in core - for beta in core - ] - ) - - if active is None: - if not core: - active = list(range(orbitals)) - else: - active = [i for i in range(orbitals) if i not in core] - - if any(i > orbitals - 1 or i < 0 for i in active): - raise ValueError( - f"Indices of active orbitals must be between 0 and {orbitals - 1}; got active = {active}" - ) - - # Indices of the matrix elements with absolute values >= cutoff - indices = np.nonzero(np.abs(matrix_elements) >= cutoff) - - # Single out the indices of active orbitals - num_indices = len(indices[0]) - quads = [ - [indices[0][i], indices[1][i], indices[2][i], indices[3][i]] - for i in range(num_indices) - if all(indices[j][i] in active for j in range(len(indices))) - ] - - # Build the FermionOperator representing V - v_op = openfermion.ops.FermionOperator("") * v_core - - # add renormalized (due to core orbitals) "one-particle" operators - if core: - for alpha in active: - for beta in active: - element = 2 * np.sum( - matrix_elements[np.array(core), alpha, beta, np.array(core)] - ) - np.sum(matrix_elements[np.array(core), alpha, np.array(core), beta]) - - # up-up term - a = 2 * active.index(alpha) - b = 2 * active.index(beta) - v_op += openfermion.ops.FermionOperator(((a, 1), (b, 0)), element) - - # down-down term - v_op += openfermion.ops.FermionOperator(((a + 1, 1), (b + 1, 0)), element) - - # add two-particle operators - for quad in quads: - alpha, beta, gamma, delta = quad - element = matrix_elements[alpha, beta, gamma, delta] - - # up-up-up-up term - a = 2 * active.index(alpha) - b = 2 * active.index(beta) - g = 2 * active.index(gamma) - d = 2 * active.index(delta) - v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) - - # up-down-down-up term - b = 2 * active.index(beta) + 1 - g = 2 * active.index(gamma) + 1 - v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) - - # down-up-up-down term - a = 2 * active.index(alpha) + 1 - b = 2 * active.index(beta) - g = 2 * active.index(gamma) - d = 2 * active.index(delta) + 1 - v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) - - # down-down-down-down term - b = 2 * active.index(beta) + 1 - g = 2 * active.index(gamma) + 1 - v_op += openfermion.ops.FermionOperator(((a, 1), (b, 1), (g, 0), (d, 0)), 0.5 * element) - - return v_op - - -def dipole_of( - symbols, - coordinates, - name="molecule", - charge=0, - mult=1, - basis="sto-3g", - package="pyscf", - core=None, - active=None, - mapping="jordan_wigner", - cutoff=1.0e-12, - outpath=".", - wires=None, -): - r"""Computes the electric dipole moment operator in the Pauli basis. - - The second quantized dipole moment operator :math:`\hat{D}` of a molecule is given by - - .. math:: - - \hat{D} = -\sum_{\alpha, \beta} \langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle - [\hat{c}_{\alpha\uparrow}^\dagger \hat{c}_{\beta\uparrow} + - \hat{c}_{\alpha\downarrow}^\dagger \hat{c}_{\beta\downarrow}] + \hat{D}_\mathrm{n}. - - In the equation above, the indices :math:`\alpha, \beta` run over the basis of Hartree-Fock - molecular orbitals and the operators :math:`\hat{c}^\dagger` and :math:`\hat{c}` are the - electron creation and annihilation operators, respectively. The matrix elements of the - position operator :math:`\hat{{\bf r}}` are computed as - - .. math:: - - \langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle = \sum_{i, j} - C_{\alpha i}^*C_{\beta j} \langle i \vert \hat{{\bf r}} \vert j \rangle, - - where :math:`\vert i \rangle` is the wave function of the atomic orbital, - :math:`C_{\alpha i}` are the coefficients defining the molecular orbitals, - and :math:`\langle i \vert \hat{{\bf r}} \vert j \rangle` - is the representation of operator :math:`\hat{{\bf r}}` in the atomic basis. - - The contribution of the nuclei to the dipole operator is given by - - .. math:: - - \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i \hat{I}, - - - where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the - nuclear coordinates of the :math:`i`-th atom of the molecule. - - Args: - symbols (list[str]): symbols of the atomic species in the molecule - coordinates (array[float]): 1D array with the atomic positions in Cartesian - coordinates. The coordinates must be given in atomic units and the size of the array - should be ``3*N`` where ``N`` is the number of atoms. - name (str): name of the molecule - charge (int): charge of the molecule - mult (int): spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1` of the - Hartree-Fock (HF) state based on the number of unpaired electrons occupying the - HF orbitals - basis (str): Atomic basis set used to represent the molecular orbitals. Basis set - availability per element can be found - `here `_ - package (str): quantum chemistry package (pyscf) used to solve the - mean field electronic structure problem - core (list): indices of core orbitals - active (list): indices of active orbitals - mapping (str): transformation (``'jordan_wigner'`` or ``'bravyi_kitaev'``) used to - map the fermionic operator to the Pauli basis - cutoff (float): Cutoff value for including the matrix elements - :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements - with absolute value less than ``cutoff`` are neglected. - outpath (str): path to the directory containing output files - wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator - to an observable measurable in a PennyLane ansatz. - For types Wires/list/tuple, each item in the iterable represents a wire label - corresponding to the qubit number equal to its index. - For type dict, only int-keyed dict (for qubit-to-wire conversion) is accepted. - If None, will use identity map (e.g. 0->0, 1->1, ...). - - Returns: - list[pennylane.Hamiltonian]: the qubit observables corresponding to the components - :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator in - atomic units. - - **Example** - - >>> symbols = ["H", "H", "H"] - >>> coordinates = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0, 1.855, 0.002, 0.0]) - >>> dipole_obs = dipole_of(symbols, coordinates, charge=1) - >>> print([(h.wires) for h in dipole_obs]) - [, , ] - - >>> dipole_obs[0] # x-component of D - ( - 0.4781123173263876 * Z(0) - + 0.4781123173263876 * Z(1) - + -0.3913638489489803 * (Y(0) @ Z(1) @ Y(2)) - + -0.3913638489489803 * (X(0) @ Z(1) @ X(2)) - + -0.3913638489489803 * (Y(1) @ Z(2) @ Y(3)) - + -0.3913638489489803 * (X(1) @ Z(2) @ X(3)) - + 0.2661114704527088 * (Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4)) - + 0.2661114704527088 * (X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4)) - + 0.2661114704527088 * (Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5)) - + 0.2661114704527088 * (X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5)) - + 0.7144779061810713 * Z(2) - + 0.7144779061810713 * Z(3) - + -0.11734958781031017 * (Y(2) @ Z(3) @ Y(4)) - + -0.11734958781031017 * (X(2) @ Z(3) @ X(4)) - + -0.11734958781031017 * (Y(3) @ Z(4) @ Y(5)) - + -0.11734958781031017 * (X(3) @ Z(4) @ X(5)) - + 0.24190977644645698 * Z(4) - + 0.24190977644645698 * Z(5) - ) - """ - openfermion, _ = _import_of() - - if mult != 1: - raise ValueError( - f"Currently, this functionality is constrained to Hartree-Fock states with spin multiplicity = 1;" - f" got multiplicity 2S+1 = {mult}" - ) - - for i in symbols: - if i not in atomic_numbers: - raise ValueError( - f"Currently, only first- or second-row elements of the periodic table are supported;" - f" got element {i}" - ) - - hf_file = qml.qchem.meanfield(symbols, coordinates, name, charge, mult, basis, package, outpath) - - hf = openfermion.MolecularData(filename=hf_file.strip()) - - # Load dipole matrix elements in the atomic basis - # pylint: disable=import-outside-toplevel - from pyscf import gto - - mol = gto.M( - atom=hf.geometry, basis=hf.basis, charge=hf.charge, spin=0.5 * (hf.multiplicity - 1) - ) - dip_ao = mol.intor_symmetric("int1e_r", comp=3).real - - # Transform dipole matrix elements to the MO basis - n_orbs = hf.n_orbitals - c_hf = hf.canonical_orbitals - - dip_mo = np.zeros((3, n_orbs, n_orbs)) - for comp in range(3): - for alpha in range(n_orbs): - for beta in range(alpha + 1): - dip_mo[comp, alpha, beta] = c_hf[:, alpha] @ dip_ao[comp] @ c_hf[:, beta] - - dip_mo[comp] += dip_mo[comp].T - np.diag(np.diag(dip_mo[comp])) - - # Compute the nuclear contribution - dip_n = np.zeros(3) - for comp in range(3): - for i, symb in enumerate(symbols): - dip_n[comp] += atomic_numbers[symb] * coordinates[3 * i + comp] - - # Build the observable - dip = [] - for i in range(3): - fermion_obs = one_particle(dip_mo[i], core=core, active=active, cutoff=cutoff) - dip.append(observable([-fermion_obs], init_term=dip_n[i], mapping=mapping, wires=wires)) - - return dip - - -def molecular_dipole( - molecule, - method="dhf", - active_electrons=None, - active_orbitals=None, - mapping="jordan_wigner", - outpath=".", - wires=None, - args=None, - cutoff=1.0e-16, -): # pylint:disable=too-many-arguments, too-many-statements - r"""Generate the dipole moment operator for a molecule in the Pauli basis. - - The dipole operator in the second-quantized form is - - .. math:: - - \hat{D} = -\sum_{pq} d_{pq} [\hat{c}_{p\uparrow}^\dagger \hat{c}_{q\uparrow} + - \hat{c}_{p\downarrow}^\dagger \hat{c}_{q\downarrow}] - - \hat{D}_\mathrm{c} + \hat{D}_\mathrm{n}, - - where the matrix elements :math:`d_{pq}` are given by the integral of the position operator - :math:`\hat{{\bf r}}` over molecular orbitals :math:`\phi` - - .. math:: - - d_{pq} = \int \phi_p^*(r) \hat{{\bf r}} \phi_q(r) dr, - - and :math:`\hat{c}^{\dagger}` and :math:`\hat{c}` are the creation and annihilation operators, - respectively. The contribution of the core orbitals and nuclei are denoted by - :math:`\hat{D}_\mathrm{c}` and :math:`\hat{D}_\mathrm{n}`, respectively, which are computed as - - .. math:: - \hat{D}_\mathrm{c} = 2 \sum_{i=1}^{N_\mathrm{core}} d_{ii} \quad \text{and} \quad - \hat{D}_\mathrm{n} = \sum_{i=1}^{N_\mathrm{atoms}} Z_i {\bf R}_i, - - where :math:`Z_i` and :math:`{\bf R}_i` denote, respectively, the atomic number and the - nuclear coordinates of the :math:`i`-th atom of the molecule. - - The fermionic dipole operator is then transformed to the qubit basis, which gives - - .. math:: - - \hat{D} = \sum_{j} c_j P_j, - - where :math:`c_j` is a numerical coefficient and :math:`P_j` is a tensor product of - single-qubit Pauli operators :math:`X, Y, Z, I`. The qubit observables corresponding - to the components :math:`\hat{D}_x`, :math:`\hat{D}_y`, and :math:`\hat{D}_z` of the - dipole operator are then computed separately. - - Args: - 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, 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 - are considered to be active. - mapping (str): Transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian. - Input values can be ``'jordan_wigner'``, ``'parity'`` or ``'bravyi_kitaev'``. - outpath (str): Path to the directory containing output files - wires (Wires, list, tuple, dict): Custom wire mapping used to convert the qubit operator to - an observable measurable in a Pennylane ansatz. - For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label - 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. - args (array[array[float]]): Initial values of the differentiable parameters - cutoff (float): Cutoff value for including the matrix elements - :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements - with absolute value less than ``cutoff`` are neglected. - - Returns: - list[pennylane.Hamiltonian]: The qubit observables corresponding to the components - :math:`\hat{D}_x`, :math:`\hat{D}_y` and :math:`\hat{D}_z` of the dipole operator. - - **Example** - - >>> symbols = ["H", "H", "H"] - >>> coordinates = np.array([[0.028, 0.054, 0.0], [0.986, 1.610, 0.0], [1.855, 0.002, 0.0]]) - >>> mol = qml.qchem.Molecule(symbols, coordinates, charge=1) - >>> dipole_obs = qml.qchem.molecular_dipole(mol, method="openfermion") - >>> dipole_obs[0] # x-component of D - ( - 0.4781123173263876 * Z(0) - + 0.4781123173263876 * Z(1) - + -0.3913638489489803 * (Y(0) @ Z(1) @ Y(2)) - + -0.3913638489489803 * (X(0) @ Z(1) @ X(2)) - + -0.3913638489489803 * (Y(1) @ Z(2) @ Y(3)) - + -0.3913638489489803 * (X(1) @ Z(2) @ X(3)) - + 0.2661114704527088 * (Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4)) - + 0.2661114704527088 * (X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4)) - + 0.2661114704527088 * (Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5)) - + 0.2661114704527088 * (X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5)) - + 0.7144779061810713 * Z(2) - + 0.7144779061810713 * Z(3) - + -0.11734958781031017 * (Y(2) @ Z(3) @ Y(4)) - + -0.11734958781031017 * (X(2) @ Z(3) @ X(4)) - + -0.11734958781031017 * (Y(3) @ Z(4) @ Y(5)) - + -0.11734958781031017 * (X(3) @ Z(4) @ X(5)) - + 0.24190977644645698 * Z(4) - + 0.24190977644645698 * Z(5) - ) - """ - - if method not in ["dhf", "openfermion"]: - raise ValueError("Only 'dhf', 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'." - ) - - symbols = molecule.symbols - coordinates = molecule.coordinates - - if qml.math.shape(coordinates)[0] == len(symbols) * 3: - geometry_dhf = qml.numpy.array(coordinates.reshape(len(symbols), 3)) - geometry_hf = coordinates - elif len(coordinates) == len(symbols): - geometry_dhf = qml.numpy.array(coordinates) - geometry_hf = coordinates.flatten() - - if molecule.mult != 1: - raise ValueError( - "Open-shell systems are not supported. Change the charge or spin multiplicity of the molecule." - ) - - core, active = qml.qchem.active_space( - molecule.n_electrons, molecule.n_orbitals, molecule.mult, active_electrons, active_orbitals - ) - - if method == "dhf": - - if args is None and isinstance(geometry_dhf, qml.numpy.tensor): - geometry_dhf.requires_grad = False - mol = qml.qchem.Molecule( - symbols, - geometry_dhf, - charge=molecule.charge, - mult=molecule.mult, - basis_name=molecule.basis_name, - load_data=molecule.load_data, - alpha=molecule.alpha, - coeff=molecule.coeff, - ) - - requires_grad = args is not None - dip = ( - qml.qchem.dipole_moment(mol, cutoff=cutoff, core=core, active=active, mapping=mapping)( - *(args or []) - ) - if requires_grad - else qml.qchem.dipole_moment( - mol, cutoff=cutoff, core=core, active=active, mapping=mapping - )() - ) - if wires: - wires_new = qml.qchem.convert._process_wires(wires) - wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels))) - dip = [qml.map_wires(op, wires_map) for op in dip] - - return dip - - dip = qml.qchem.dipole_of( - symbols, - geometry_hf, - molecule.name, - molecule.charge, - molecule.mult, - molecule.basis_name, - package="pyscf", - core=core, - active=active, - mapping=mapping, - cutoff=cutoff, - outpath=outpath, - wires=wires, - ) - - return dip - - -def meanfield( - symbols, - coordinates, - name="molecule", - charge=0, - mult=1, - basis="sto-3g", - package="pyscf", - outpath=".", -): # pylint: disable=too-many-arguments - r"""Generates a file from which the mean field electronic structure - of the molecule can be retrieved. - - This function uses OpenFermion-PySCF plugins to - perform the Hartree-Fock (HF) calculation for the polyatomic system using the quantum - chemistry packages ``PySCF``. The mean field electronic - structure is saved in an hdf5-formatted file. - - The charge of the molecule can be given to simulate cationic/anionic systems. - Also, the spin multiplicity can be input to determine the number of unpaired electrons - occupying the HF orbitals as illustrated in the figure below. - - | - - .. figure:: ../../_static/qchem/hf_references.png - :align: center - :width: 50% - - | - - Args: - symbols (list[str]): symbols of the atomic species in the molecule - coordinates (array[float]): 1D array with the atomic positions in Cartesian - coordinates. The coordinates must be given in atomic units and the size of the array - should be ``3*N`` where ``N`` is the number of atoms. - name (str): molecule label - charge (int): net charge of the system - mult (int): Spin multiplicity :math:`\mathrm{mult}=N_\mathrm{unpaired} + 1` for - :math:`N_\mathrm{unpaired}` unpaired electrons occupying the HF orbitals. - Possible values for ``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 HF orbitals. Basis set - availability per element can be found - `here `_ - package (str): Quantum chemistry package used to solve the Hartree-Fock equations. - outpath (str): path to output directory - - Returns: - str: absolute path to the file containing the mean field electronic structure - - **Example** - - >>> symbols, coordinates = (['H', 'H'], np.array([0., 0., -0.66140414, 0., 0., 0.66140414])) - >>> meanfield(symbols, coordinates, name="h2") - ./h2_pyscf_sto-3g - """ - openfermion, openfermionpyscf = _import_of() - - if coordinates.size != 3 * len(symbols): - raise ValueError( - f"The size of the array 'coordinates' has to be 3*len(symbols) = {3 * len(symbols)};" - f" got 'coordinates.size' = {coordinates.size}" - ) - - package = package.strip().lower() - - if package not in "pyscf": - error_message = ( - f"Integration with quantum chemistry package '{package}' is not available. \n Please set" - f" 'package' to 'pyscf'." - ) - raise TypeError(error_message) - - filename = name + "_" + package.lower() + "_" + basis.strip() - path_to_file = os.path.join(outpath.strip(), filename) - - geometry = [ - [symbol, tuple(np.array(coordinates)[3 * i : 3 * i + 3] * bohr_angs)] - for i, symbol in enumerate(symbols) - ] - - molecule = openfermion.MolecularData(geometry, basis, mult, charge, filename=path_to_file) - - if package == "pyscf": - # pylint: disable=import-outside-toplevel - from openfermionpyscf import run_pyscf - - run_pyscf(molecule, run_scf=1, verbose=0) - - return path_to_file - - -def decompose(hf_file, mapping="jordan_wigner", core=None, active=None): - r"""Decomposes the molecular Hamiltonian into a linear combination of Pauli operators using - OpenFermion tools. - - This function uses OpenFermion functions to build the second-quantized electronic Hamiltonian - of the molecule and map it to the Pauli basis using the Jordan-Wigner or Bravyi-Kitaev - transformation. - - Args: - hf_file (str): absolute path to the hdf5-formatted file with the - Hartree-Fock electronic structure - mapping (str): Specifies the transformation to map the fermionic Hamiltonian to the - Pauli basis. Input values can be ``'jordan_wigner'`` or ``'bravyi_kitaev'``. - core (list): indices of core orbitals, i.e., the orbitals that are - not correlated in the many-body wave function - active (list): indices of active orbitals, i.e., the orbitals used to - build the correlated many-body wave function - - Returns: - QubitOperator: an instance of OpenFermion's ``QubitOperator`` - - **Example** - - >>> decompose('./pyscf/sto-3g/h2', mapping='bravyi_kitaev') - (-0.04207897696293986+0j) [] + (0.04475014401986122+0j) [X0 Z1 X2] + - (0.04475014401986122+0j) [X0 Z1 X2 Z3] +(0.04475014401986122+0j) [Y0 Z1 Y2] + - (0.04475014401986122+0j) [Y0 Z1 Y2 Z3] +(0.17771287459806262+0j) [Z0] + - (0.17771287459806265+0j) [Z0 Z1] +(0.1676831945625423+0j) [Z0 Z1 Z2] + - (0.1676831945625423+0j) [Z0 Z1 Z2 Z3] +(0.12293305054268105+0j) [Z0 Z2] + - (0.12293305054268105+0j) [Z0 Z2 Z3] +(0.1705973832722409+0j) [Z1] + - (-0.2427428049645989+0j) [Z1 Z2 Z3] +(0.1762764080276107+0j) [Z1 Z3] + - (-0.2427428049645989+0j) [Z2] - """ - openfermion, _ = _import_of() - - # loading HF data from the hdf5 file - molecule = openfermion.MolecularData(filename=hf_file.strip()) - - # getting the terms entering the second-quantized Hamiltonian - terms_molecular_hamiltonian = molecule.get_molecular_hamiltonian( - occupied_indices=core, active_indices=active - ) - - # generating the fermionic Hamiltonian - fermionic_hamiltonian = openfermion.transforms.get_fermion_operator(terms_molecular_hamiltonian) - - mapping = mapping.strip().lower() - - # 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(*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. - - The net charge of the molecule can be given to simulate cationic/anionic systems. Also, the - spin multiplicity can be input to determine the number of unpaired electrons occupying the HF - orbitals as illustrated in the left panel of the figure below. - - The basis of Gaussian-type *atomic* orbitals used to represent the *molecular* orbitals can be - specified to go beyond the minimum basis approximation. - - An active space can be defined for a given number of *active electrons* occupying a reduced set - of *active orbitals* as sketched in the right panel of the figure below. - - | - - .. figure:: ../../_static/qchem/fig_mult_active_space.png - :align: center - :width: 90% - - | - - Args: - 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 - 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 - are considered to be active. - mapping (str): transformation used to map the fermionic Hamiltonian to the qubit Hamiltonian - outpath (str): path to the directory containing output files - wires (Wires, list, tuple, dict): Custom wire mapping for connecting to Pennylane ansatz. - For types ``Wires``/``list``/``tuple``, each item in the iterable represents a wire label - 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. - args (array[array[float]]): initial values of the differentiable parameters - convert_tol (float): Tolerance in `machine epsilon `_ - 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 - 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 = ['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) - (-0.04207897647782188) [I0] - + (0.17771287465139934) [Z0] - + (0.1777128746513993) [Z1] - + (-0.24274280513140484) [Z2] - + (-0.24274280513140484) [Z3] - + (0.17059738328801055) [Z0 Z1] - + (0.04475014401535161) [Y0 X1 X2 Y3] - + (-0.04475014401535161) [Y0 Y1 X2 X3] - + (-0.04475014401535161) [X0 X1 Y2 Y3] - + (0.04475014401535161) [X0 Y1 Y2 X3] - + (0.12293305056183801) [Z0 Z2] - + (0.1676831945771896) [Z0 Z3] - + (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 - elif len(coordinates) == len(symbols): - 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 in ("dhf", "pyscf"): - n_electrons = sum([atomic_numbers[s] for s in symbols]) - charge - - if n_electrons % 2 == 1 or mult != 1: - raise ValueError( - "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( - symbols, - geometry_dhf, - charge=charge, - mult=mult, - basis_name=basis, - load_data=load_data, - alpha=alpha, - coeff=coeff, - ) - core, active = qml.qchem.active_space( - mol.n_electrons, mol.n_orbitals, mult, active_electrons, active_orbitals - ) - - requires_grad = args is not None - h = ( - 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, mapping=mapping)() - ) - - 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: - coeffs = qml.numpy.real(h.coeffs, requires_grad=requires_grad) - h = qml.Hamiltonian(coeffs, h.ops) - - if wires: - h = qml.map_wires(h, wires_map) - return h, 2 * len(active) - - 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(): - 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: - 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) - - openfermion, _ = _import_of() - - hf_file = meanfield(symbols, geometry_hf, name, charge, mult, basis, "pyscf", outpath) - - molecule = openfermion.MolecularData(filename=hf_file) - - core, active = qml.qchem.active_space( - molecule.n_electrons, molecule.n_orbitals, mult, active_electrons, active_orbitals - ) - - h_of, qubits = (decompose(hf_file, mapping, core, active), 2 * len(active)) - - h_pl = qml.qchem.convert.import_operator(h_of, wires=wires, tol=convert_tol) - - 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 From fe779b6115c819d65a6608952fb7dca77b231c10 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Thu, 6 Jun 2024 13:01:16 -0400 Subject: [PATCH 27/28] Fixed CI --- pennylane/qchem/openfermion_pyscf.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pennylane/qchem/openfermion_pyscf.py b/pennylane/qchem/openfermion_pyscf.py index 76e878b82d9..bc396b37fbe 100644 --- a/pennylane/qchem/openfermion_pyscf.py +++ b/pennylane/qchem/openfermion_pyscf.py @@ -653,6 +653,7 @@ def dipole_of( return dip + def molecular_dipole( molecule, method="dhf", @@ -937,6 +938,7 @@ def decompose(hf_file, mapping="jordan_wigner", core=None, active=None): This function uses OpenFermion functions to build the second-quantized electronic Hamiltonian of the molecule and map it to the Pauli basis using the Jordan-Wigner, Parity or Bravyi-Kitaev + transformation. Args: hf_file (str): absolute path to the hdf5-formatted file with the From 5b364d259e397aa83b6ec2f6c9174be312ee97f0 Mon Sep 17 00:00:00 2001 From: ddhawan11 Date: Fri, 7 Jun 2024 06:03:05 -0400 Subject: [PATCH 28/28] change docstring to include parity --- pennylane/qchem/openfermion_pyscf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pennylane/qchem/openfermion_pyscf.py b/pennylane/qchem/openfermion_pyscf.py index bc396b37fbe..91d2faf8bc3 100644 --- a/pennylane/qchem/openfermion_pyscf.py +++ b/pennylane/qchem/openfermion_pyscf.py @@ -98,7 +98,7 @@ def observable(fermion_ops, init_term=0, mapping="jordan_wigner", wires=None): - The function uses tools of `OpenFermion `_ to map the resulting fermionic Hamiltonian to the basis of Pauli matrices via the - Jordan-Wigner or Bravyi-Kitaev transformation. Finally, the qubit operator is converted + Jordan-Wigner, Parity or Bravyi-Kitaev transformation. Finally, the qubit operator is converted to a PennyLane observable by the function :func:`~.convert_observable`. Args: @@ -551,7 +551,7 @@ def dipole_of( mean field electronic structure problem core (list): indices of core orbitals active (list): indices of active orbitals - mapping (str): transformation (``'jordan_wigner'`` or ``'bravyi_kitaev'``) used to + mapping (str): transformation (``'jordan_wigner'``, ``'parity'``, or ``'bravyi_kitaev'``) used to map the fermionic operator to the Pauli basis cutoff (float): Cutoff value for including the matrix elements :math:`\langle \alpha \vert \hat{{\bf r}} \vert \beta \rangle`. The matrix elements