diff --git a/circuit_knitting/forging/entanglement_forging_ground_state_solver.py b/circuit_knitting/forging/entanglement_forging_ground_state_solver.py index 872bf2728..9baffe5f3 100644 --- a/circuit_knitting/forging/entanglement_forging_ground_state_solver.py +++ b/circuit_knitting/forging/entanglement_forging_ground_state_solver.py @@ -129,6 +129,7 @@ def __init__( backend_names: str | list[str] | None = None, options: Options | list[Options] | None = None, mo_coeff: np.ndarray | None = None, + hf_energy: float | None = None, ): """ Assign the necessary class variables and initialize any defaults. @@ -143,6 +144,9 @@ def __init__( backend_names: Backend name or list of backend names to use during parallel computation options: Options or list of options to be applied to the backends mo_coeff: Coefficients for converting an input problem to MO basis + hf_energy: The Hartree-Fock energy to use at each iteration. If this is left unset, + the energy will be calculated using quantum experiments. See :ref:`Fixing the Hartree-Fock bitstring` + for more information. Returns: None @@ -157,6 +161,7 @@ def __init__( self._orbitals_to_reduce = orbitals_to_reduce self.backend_names = backend_names # type: ignore self.options = options + self.hf_energy = hf_energy self._mo_coeff = mo_coeff self._optimizer: Optimizer | MINIMIZER = optimizer or SPSA() @@ -246,6 +251,16 @@ def mo_coeff(self, mo_coeff: np.ndarray | None) -> None: """Set the coefficients for converting integrals to the MO basis.""" self._mo_coeff = mo_coeff + @property + def hf_energy(self) -> float | None: + """Return the Hartree-Fock energy.""" + return self._hf_energy + + @hf_energy.setter + def hf_energy(self, hf_energy: float | None) -> None: + """Set the Hartree-Fock energy.""" + self._hf_energy = hf_energy + def solve( self, problem: ElectronicStructureProblem, @@ -284,18 +299,24 @@ def solve( hamiltonian_terms = self.get_qubit_operators(problem) ef_operator = convert_cholesky_operator(hamiltonian_terms, self._ansatz) - # Set the knitter class field + # Shift the hf value so it can be entered into the Schmidt matrix + fixed_hf_value = None + if self.hf_energy is not None: + fixed_hf_value = self.hf_energy - self._energy_shift if self._service is not None: backend_names = self._backend_names or ["ibmq_qasm_simulator"] self._knitter = EntanglementForgingKnitter( self._ansatz, + fixed_hf_value=fixed_hf_value, service=self._service, backend_names=backend_names, options=self._options, ) else: self._knitter = EntanglementForgingKnitter( - self._ansatz, options=self._options + self._ansatz, + fixed_hf_value=fixed_hf_value, + options=self._options, ) self._history = EntanglementForgingHistory() self._eval_count = 0 diff --git a/circuit_knitting/forging/entanglement_forging_knitter.py b/circuit_knitting/forging/entanglement_forging_knitter.py index 894b5a48e..612923b31 100644 --- a/circuit_knitting/forging/entanglement_forging_knitter.py +++ b/circuit_knitting/forging/entanglement_forging_knitter.py @@ -13,6 +13,7 @@ from __future__ import annotations +import copy from typing import Sequence, Any from concurrent.futures import ThreadPoolExecutor @@ -38,6 +39,7 @@ class EntanglementForgingKnitter: def __init__( self, ansatz: EntanglementForgingAnsatz, + fixed_hf_value: float | None = None, service: QiskitRuntimeService | None = None, backend_names: str | list[str] | None = None, options: Options | list[Options] | None = None, @@ -48,6 +50,9 @@ def __init__( Args: ansatz: The container for the circuit structure and bitstrings to be used (and generate the stateprep circuits) + fixed_hf_value: This value will be used instead of calculating the Hartree-Fock energy. + This value should represent the Hartree-Fock energy shifted by the nuclear repulsion energy. + The HF energy needs to be shifted in order to be used in the Schmidt matrix. service: The service used to spawn Qiskit primitives and runtime jobs backend_names: Names of the backends to use for calculating expectation values options: Options to use with the backends @@ -59,6 +64,7 @@ def __init__( self._session_ids: list[str | None] | None = None self._backend_names: list[str] | None = None self._options: list[Options] | None = None + self.fixed_hf_value = fixed_hf_value # Call constructors self.backend_names = backend_names # type: ignore self.options = options @@ -141,6 +147,16 @@ def service(self, service: QiskitRuntimeService | None) -> None: """Change the service class field.""" self._service = service.active_account() if service is not None else service + @property + def fixed_hf_value(self) -> float | None: + """Return the shifted Hartree-Fock energy.""" + return self._fixed_hf_value + + @fixed_hf_value.setter + def fixed_hf_value(self, fixed_hf_value: float | None) -> None: + """Set the shifted Hartree-Fock energy to bet used in the Schmidt matrix.""" + self._fixed_hf_value = fixed_hf_value + def __call__( self, ansatz_parameters: Sequence[float], @@ -191,6 +207,8 @@ def __call__( tensor_ansatze_u = [ prep_circ.compose(circuit_u) for prep_circ in self._tensor_circuits_u ] + if self.fixed_hf_value is not None: + tensor_ansatze_u = tensor_ansatze_u[1:] superposition_ansatze_u = [ prep_circ.compose(circuit_u) for prep_circ in self._superposition_circuits_u ] @@ -201,6 +219,8 @@ def __call__( tensor_ansatze_v = [ prep_circ.compose(circuit_u) for prep_circ in self._tensor_circuits_v ] + if self.fixed_hf_value is not None: + tensor_ansatze_v = tensor_ansatze_v[1:] superposition_ansatze_v = [ prep_circ.compose(circuit_u) for prep_circ in self._superposition_circuits_v @@ -282,10 +302,29 @@ def __call__( ) self._session_ids[i] = job_id + # Pack NaNs into the expval list in place of the HF bitstring results to + # ensure Schmidt matrix is the correct shape + if self._fixed_hf_value is not None: + num_paulis = len(forged_operator.tensor_paulis) + nans_u = np.empty(num_paulis) + nans_u[:] = np.nan + if not self._ansatz.bitstrings_are_symmetric: + # Should be equal number of expvals for each subsystem + assert len(tensor_expvals) % 2 == 0 + num_tensor_terms = int(len(tensor_expvals) / 2) + nans_v = copy.deepcopy(nans_u) + tensor_expvals.insert(num_tensor_terms, nans_v) + tensor_expvals.insert(0, nans_u) + # Compute the Schmidt matrix h_schmidt = self._compute_h_schmidt( forged_operator, np.array(tensor_expvals), np.array(superposition_expvals) ) + + # Hard-code the shifted Hartree-Fock energy, if desired + if self._fixed_hf_value is not None: + h_schmidt[0, 0] = self._fixed_hf_value + evals, evecs = np.linalg.eigh(h_schmidt) schmidt_coeffs = evecs[:, 0] energy = evals[0] @@ -569,10 +608,10 @@ def _estimate_expvals( tensor_paulis: list[Pauli], superposition_ansatze: list[QuantumCircuit], superposition_paulis: list[Pauli], - service_args: dict[str, Any] | None = None, - backend_name: str | None = None, - options: Options | None = None, - session_id: str | None = None, + service_args: dict[str, Any] | None, + backend_name: str | None, + options: Options | None, + session_id: str | None, ) -> tuple[list[np.ndarray], list[np.ndarray], str | None]: """Run quantum circuits to generate the expectation values. diff --git a/docs/entanglement_forging/explanation/figs/fixed_hf.png b/docs/entanglement_forging/explanation/figs/fixed_hf.png new file mode 100644 index 000000000..7df1d0b64 Binary files /dev/null and b/docs/entanglement_forging/explanation/figs/fixed_hf.png differ diff --git a/docs/entanglement_forging/explanation/index.rst b/docs/entanglement_forging/explanation/index.rst index a89773bc0..9f77bbe4e 100644 --- a/docs/entanglement_forging/explanation/index.rst +++ b/docs/entanglement_forging/explanation/index.rst @@ -177,6 +177,44 @@ that do not participate in electronic excitations (i.e. core orbitals or those that lie out of symmetry) by removing the bits that correspond to them. +.. _Fixing the Hartree-Fock bitstring: + +Fixing the Hartree-Fock bitstring +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In some cases, it is possible to increase the accuracy of simulations and speed up +the execution by bypassing the experiments associated with the first bitstring and +replacing those results with the Hartree-Fock energy value. + +.. code-block:: python + :caption: Fixing the HF energy at each iteration + from qiskit_nature.second_q.problems import ElectronicStructureProblem + from circuit_knitting.forging import EntanglementForgingGroundStateSolver + + problem = ElectronicStructureProblem(...) + hf_energy = ... + + solver = EntanglementForgingGroundStateSolver( + ansatz=ansatz, + hf_energy=hf_energy + ) + + result = solver.solve(problem) + +This setting requires an ansatz that leaves the Hartree-Fock (HF) state +unchanged with respect to the optimization parameters. As a rule of thumb, +this can be achieved by restricting entanglement between the qubits representing +occupied orbitals (bits = 1) in the HF state and the qubits representing +unoccupied orbitals (bits = 0) in the HF state. + +For example, this figure from [1] shows the A, B, and C qubits entangled with +the hop gates, D & E qubits entangled with hop gates, while the partition between +(A,B,C) and (D,E) are only entangled with a CZ gate. + +.. figure:: figs/fixed_hf.png + :width: 250 + :alt: Fixing the first bitstring to the HF value + .. _Ansatz design: Designing the ansatz used in Entanglement Forging diff --git a/releasenotes/notes/fixed-hartree-fock-447c33a956ffea84.yaml b/releasenotes/notes/fixed-hartree-fock-447c33a956ffea84.yaml new file mode 100644 index 000000000..c6b56aece --- /dev/null +++ b/releasenotes/notes/fixed-hartree-fock-447c33a956ffea84.yaml @@ -0,0 +1,4 @@ +--- +features: + - | + Users may now bypass experiments associated with the Hartree-Fock bitstring and replace their results with a specified Hartree-Fock value using the ``hf_energy`` class field in :class:`~circuit_knitting.forging.EntanglementForgingGroundStateSolver`. Refer to the :ref:`explanatory material ` for more information. diff --git a/test/forging/test_data/H2O_0.90_one_body.npy b/test/forging/test_data/H2O_0.90_one_body.npy new file mode 100644 index 000000000..3d48f83eb Binary files /dev/null and b/test/forging/test_data/H2O_0.90_one_body.npy differ diff --git a/test/forging/test_data/H2O_0.90_two_body.npy b/test/forging/test_data/H2O_0.90_two_body.npy new file mode 100644 index 000000000..48c28c840 Binary files /dev/null and b/test/forging/test_data/H2O_0.90_two_body.npy differ diff --git a/test/forging/test_entanglement_forging_ground_state_solver.py b/test/forging/test_entanglement_forging_ground_state_solver.py index dd6060fd5..a9b6f3eed 100644 --- a/test/forging/test_entanglement_forging_ground_state_solver.py +++ b/test/forging/test_entanglement_forging_ground_state_solver.py @@ -11,18 +11,25 @@ """Tests for EntanglementForgingGroundStateSolver module.""" +import os import unittest import importlib.util -import numpy as np +import pytest +import numpy as np from qiskit.algorithms.optimizers import SPSA +from qiskit.circuit import QuantumCircuit, Parameter from qiskit.circuit.library import TwoLocal +from qiskit.algorithms.optimizers import COBYLA from qiskit_nature.units import DistanceUnit from qiskit_nature.second_q.drivers import PySCFDriver from qiskit_nature.second_q.problems import ( ElectronicBasis, + ElectronicStructureProblem, ) +from qiskit_nature.second_q.transformers import ActiveSpaceTransformer from qiskit_nature.second_q.formats import get_ao_to_mo_from_qcschema +from qiskit_nature.second_q.hamiltonians import ElectronicEnergy from circuit_knitting.forging import ( EntanglementForgingAnsatz, @@ -35,7 +42,54 @@ class TestEntanglementForgingGroundStateSolver(unittest.TestCase): def setUp(self): # Hard-code some ansatz params/lambdas - self.optimizer = SPSA(maxiter=0) + self.hcore = np.load( + os.path.join( + os.path.dirname(__file__), "test_data", "H2O_0.90_one_body.npy" + ), + ) + self.eri = np.load( + os.path.join( + os.path.dirname(__file__), "test_data", "H2O_0.90_two_body.npy" + ), + ) + orb_act = list(range(0, 5)) + num_alpha = num_beta = 3 + hamiltonian = ElectronicEnergy.from_raw_integrals(self.hcore, self.eri) + hamiltonian.nuclear_repulsion_energy = -61.57756706745154 + problem = ElectronicStructureProblem(hamiltonian) + problem.basis = ElectronicBasis.MO + problem.num_particles = (num_alpha, num_beta) + transformer = ActiveSpaceTransformer( + num_electrons=6, num_spatial_orbitals=len(orb_act), active_orbitals=orb_act + ) + self.problem_reduced = transformer.transform(problem) + + theta = Parameter("θ") + + hop_gate = QuantumCircuit(2, name="hop_gate") + hop_gate.h(0) + hop_gate.cx(1, 0) + hop_gate.cx(0, 1) + hop_gate.ry(-theta, 0) + hop_gate.ry(-theta, 1) + hop_gate.cx(0, 1) + hop_gate.h(0) + + theta_1, theta_2, theta_3, theta_4 = ( + Parameter("θ1"), + Parameter("θ2"), + Parameter("θ3"), + Parameter("θ4"), + ) + + circuit = QuantumCircuit(5) + circuit.append(hop_gate.to_gate({theta: theta_1}), [0, 1]) + circuit.append(hop_gate.to_gate({theta: theta_2}), [3, 4]) + circuit.append(hop_gate.to_gate({theta: 0}), [1, 4]) + circuit.append(hop_gate.to_gate({theta: theta_3}), [0, 2]) + circuit.append(hop_gate.to_gate({theta: theta_4}), [3, 4]) + + self.circuit = circuit @unittest.skipIf(not pyscf_available, "pyscf is not installed") def test_entanglement_forging_vqe_hydrogen(self): @@ -62,7 +116,7 @@ def test_entanglement_forging_vqe_hydrogen(self): # Set up the entanglement forging vqe object solver = EntanglementForgingGroundStateSolver( ansatz=ansatz, - optimizer=self.optimizer, + optimizer=SPSA(maxiter=0), initial_point=[0.0, np.pi / 2], mo_coeff=mo_coeff, ) @@ -73,3 +127,77 @@ def test_entanglement_forging_vqe_hydrogen(self): # Ensure ground state energy output is within tolerance self.assertAlmostEqual(ground_state_energy, -1.121936544469326) + + @unittest.skipIf(not pyscf_available, "pyscf is not installed") + def test_fixed_hf_h2o(self): + """ + Test for fixing the HF value in computing the energy of a H2O molecule. + + Hard-coded values were generated using PySCF. + """ + # Set up the ElectronicStructureProblem + HF = -14.09259461609392 + + bitstrings = [(1, 1, 1, 0, 0), (0, 1, 1, 0, 1), (1, 1, 0, 1, 0)] + ansatz = EntanglementForgingAnsatz( + circuit_u=self.circuit, + bitstrings_u=bitstrings, + ) + + COBYLA(maxiter=0) + initial_point = [-0.83604922, -0.87326138, -0.93964018, 0.55224467] + solver = EntanglementForgingGroundStateSolver( + ansatz=ansatz, + optimizer=COBYLA(maxiter=0), + hf_energy=HF, + initial_point=initial_point, + ) + result = solver.solve(self.problem_reduced) + + assert np.allclose( + result.groundstate[0], [-0.00252657, 0.99945784, -0.03282741], atol=1e-8 + ) + + @pytest.mark.slow + @unittest.skipIf(not pyscf_available, "pyscf is not installed") + def test_fixed_hf_h2o_asymmetric(self): + """ + Test for fixing the HF value in two separate subsystems. + + Hard-coded values were generated using PySCF. + """ + # The hard-coded HF value based off pyscf outputs, given the + # molecular setup represented by the bitstrings and the problem. + HF = -14.09259461609392 + bitstrings_u = [ + (1, 1, 1, 0, 0), + (0, 1, 1, 0, 1), + (1, 1, 0, 1, 0), + (1, 1, 0, 1, 0), + ] + bitstrings_v = [ + (1, 1, 1, 0, 0), + (0, 1, 1, 0, 1), + (1, 1, 0, 1, 0), + (1, 1, 0, 0, 1), + ] + ansatz = EntanglementForgingAnsatz( + circuit_u=self.circuit, + bitstrings_u=bitstrings_u, + bitstrings_v=bitstrings_v, + ) + # The initial point was found by doing an ansatz parameter optimization using + # the reference forging implementation + initial_point = [-0.83604922, -0.87326138, -0.93964018, 0.55224467] + solver = EntanglementForgingGroundStateSolver( + ansatz=ansatz, + optimizer=COBYLA(maxiter=0), + hf_energy=HF, + initial_point=initial_point, + ) + result = solver.solve(self.problem_reduced) + + # Make sure our masks didn't interfere with the schmidt matrix + # in an unintended way. + assert not np.isnan(result.groundstate[0]).any() + assert not np.isnan(result.groundenergy)