Fix expectation with measurements code

doc/source/api-reference/ansatz.rst

@@ -19,3 +19,9 @@ Unitary Coupled Cluster
 -----------------------
 
 .. autofunction:: qibochem.ansatz.ucc.ucc_circuit
+
+
+Basis rotation
+--------------
+
+.. autofunction:: qibochem.ansatz.basis_rotation.br_circuit

doc/source/api-reference/measurement.rst

@@ -1,2 +1,6 @@
 Measurement
 ===========
+
+This section covers the API reference for obtaining the expectation value of some (molecular) Hamiltonian
+
+.. autofunction:: qibochem.measurement.expectation.expectation

pyproject.toml

@@ -6,7 +6,7 @@ authors = ["The Qibo team"]
 license = "Apache License 2.0"
 readme = "README.md"
 repository = "https://github.com/qiboteam/qibochem/"
-documentation = ""
+documentation = "https://qibo.science/qibochem/latest/"
 keywords = []
 classifiers = [
   "Programming Language :: Python :: 3",

src/qibochem/ansatz/basis_rotation.py

@@ -118,15 +118,20 @@ def givens_rotation_gate(n_qubits, orb1, orb2, theta):
 
 
 def br_circuit(n_qubits, parameters, n_occ):
-    r"""
-    Google's basis rotation circuit between the occupied/virtual orbitals. Forms the \exp(\kappa) matrix, decomposes
-        it into Givens rotations, and sets the circuit parameters based on the Givens rotation decomposition
+    """
+    Google's basis rotation circuit, applied between the occupied/virtual orbitals. Forms the exp(kappa) matrix, decomposes
+    it into Givens rotations, and sets the circuit parameters based on the Givens rotation decomposition. Note: Supposed
+    to be used with the JW fermion-to-qubit mapping
 
     Args:
         n_qubits: Number of qubits in the quantum circuit
-        parameters: Rotation parameters for \exp(\kappa)
+        parameters: Rotation parameters for exp(kappa); Must have (n_occ * n_virt) parameters
         n_occ: Number of occupied orbitals
+
+    Returns:
+        Qibo ``Circuit`` corresponding to the basis rotation ansatz between the occupied and virtual orbitals
     """
+    assert len(parameters) == (n_occ * (n_qubits - n_occ)), "Need len(parameters) == (n_occ * n_virt)"
     # Unitary rotation matrix \exp(\kappa)
     exp_k = unitary_rot_matrix(parameters, range(n_occ), range(n_occ, n_qubits))
     # List of Givens rotation parameters

src/qibochem/ansatz/hardware_efficient.py

@@ -2,8 +2,9 @@
 
 
 def hea(n_layers, n_qubits, parameter_gates=["RY", "RZ"], coupling_gates="CZ"):
-    """Builds the generalized hardware-efficient ansatz, i.e. the choice of rotation and entangling gates can be
-    manually defined
+    """
+    Builds the generalized hardware-efficient ansatz, in which the rotation and entangling gates used can be
+    chosen by the user
 
     Args:
         n_layers: Number of layers of rotation and entangling gates

src/qibochem/measurement/expectation.py

@@ -1,118 +1,57 @@
-# from qibo import gates
-# from qibo.models import Circuit
-
-import qibo
-from qibo.hamiltonians import SymbolicHamiltonian
-
-from qibochem.measurement.basis_rotate import measure_rotate_basis
-
-
-def pauli_expectation_shots(qc, nshots):
-    """
-    Calculates expectation value of circuit for pauli string from shots
-
-    Args:
-        qc: quantum circuit with appropriate rotations and measure gates
-        nshots: number of shots
-
-    Returns:
-        expectation: expectation value
-    """
-    result = qc.execute(nshots=nshots)
-    meas = result.frequencies(binary=True)
-
-    bitcount = 0
-    for bitstring, count in meas.items():
-        if bitstring.count("1") % 2 == 0:
-            bitcount += count
-        else:
-            bitcount -= count
-    expectation = bitcount / nshots
-    return expectation
-
-
-def circuit_expectation_shots(qc, of_qubham, nshots=1000):
-    """
-    Calculates expectation value of qubit hamiltonian w.r.t. circuit ansatz using shots
-
-    Args:
-        qc: Quantum circuit with ansatz
-        of_qubham: Molecular Hamiltonian given as an OpenFermion QubitOperator
-
-    Returns:
-        Expectation value of of_qubham w.r.t. circuit ansatz
-    """
-    # nterms = len(of_qubham.terms)
-    nqubits = qc.nqubits
-    # print(nterms)
-    coeffs = []
-    strings = []
-
-    expectation = 0
-    for qubop in of_qubham.terms:
-        coeffs.append(of_qubham.terms[qubop])
-        strings.append([qubop])
-
-    istring = 0
-    for pauli_op in strings:
-        if len(pauli_op[0]) == 0:  # no pauli obs to measure,
-            expectation += coeffs[istring]
-            # print(coeffs[istring])
-        else:
-            # add rotation gates to rotate pauli observable to computational basis
-            # i.e. X to Z, Y to Z
-            meas_qc = measure_rotate_basis(pauli_op[0], nqubits)
-            full_qc = qc + meas_qc
-            # print(full_qc.draw())
-            pauli_op_exp = pauli_expectation_shots(full_qc, nshots)
-            expectation += coeffs[istring] * pauli_op_exp
-            ## print(pauli_op, pauli_op_exp, coeffs[istring])
-        istring += 1
-
-    return expectation
-
-
-def expectation(
-    circuit: qibo.models.Circuit, hamiltonian: SymbolicHamiltonian, from_samples=False, n_shots=1000
-) -> float:
-    """
-    Calculate expectation value of Hamiltonian using either the state vector from running a
-        circuit, or the frequencies of the resultant binary string results
-
-    Args:
-        circuit (qibo.models.Circuit): Quantum circuit ansatz
-        hamiltonian (SymbolicHamiltonian): Molecular Hamiltonian
-        from_samples: Whether the expectation value calculation uses samples or the simulated
-            state vector. Default: False, state vector simulation
-        n_shots: Number of times the circuit is run for the from_samples=True case
-
-    Returns:
-        Hamiltonian expectation value (float)
-    """
-    if from_samples:
-        raise NotImplementedError("expectation function only works with state vector")
-    # TODO: Rough code for expectation_from_samples if issue resolved
-    # Yet to test!!!!!
-    #
-    # from functools import reduce
-    # total = 0.0
-    # Iterate over each term in the Hamiltonian
-    # for term in hamiltonian.terms:
-    #     # Get the basis rotation gates and target qubits from the Hamiltonian term
-    #     qubits = [factor.target_qubit for factor in term.factors]
-    #     basis = [type(factor.gate) for factor in term.factors]
-    #     # Run a copy of the initial circuit to get the output frequencies
-    #     _circuit = circuit.copy()
-    #     _circuit.add(gates.M(*qubits, basis=basis))
-    #     result = _circuit(nshots=n_shots)
-    #     frequencies = result.frequencies(binary=True)
-    #     # Only works for Z terms, raises an error if ham_term has X/Y terms
-    #     total += SymbolicHamiltonian(
-    #                  reduce(lambda x, y: x*y, term.factors, 1)
-    #              ).expectation_from_samples(frequencies, qubit_map=qubits)
-    # return total
-
-    # Expectation value from state vector simulation
-    result = circuit(nshots=1)
-    state_ket = result.state()
-    return hamiltonian.expectation(state_ket)
+import qibo
+from qibo import gates
+from qibo.hamiltonians import SymbolicHamiltonian
+from qibo.symbols import Z
+
+
+def expectation(
+    circuit: qibo.models.Circuit, hamiltonian: SymbolicHamiltonian, from_samples=False, n_shots=1000
+) -> float:
+    """
+    Calculate expectation value of some Hamiltonian using either the state vector or sample measurements from running a
+    quantum circuit
+
+    Args:
+        circuit (qibo.models.Circuit): Quantum circuit ansatz
+        hamiltonian (SymbolicHamiltonian): Molecular Hamiltonian
+        from_samples (Boolean): Whether the expectation value calculation uses samples or the simulated
+            state vector. Default: ``False``; Results are from a state vector simulation
+        n_shots (int): Number of times the circuit is run if ``from_samples=True``. Default: ``1000``
+
+    Returns:
+        Hamiltonian expectation value (float)
+    """
+    if from_samples:
+        from functools import reduce
+
+        total = 0.0
+        # Iterate over each term in the Hamiltonian
+        for term in hamiltonian.terms:
+            # Get the target qubits and basis rotation gates from the Hamiltonian term
+            qubits = [int(factor.target_qubit) for factor in term.factors]
+            basis = [type(factor.gate) for factor in term.factors]
+            # Run a copy of the original circuit to get the output frequencies
+            _circuit = circuit.copy()
+            _circuit.add(gates.M(*qubits, basis=basis))
+            result = _circuit(nshots=n_shots)
+            frequencies = result.frequencies(binary=True)
+            # Only works for Z terms, raises an error if ham_term has X/Y terms
+            # total += SymbolicHamiltonian(
+            #              reduce(lambda x, y: x*y, term.factors, 1)
+            #          ).expectation_from_samples(frequencies, qubit_map=qubits)
+            # Workaround code to handle X/Y terms in the Hamiltonian:
+            # Get each Pauli string e.g. X0Y1
+            pauli = [factor.name for factor in term.factors]
+            # Replace each X and Y symbol with Z; then include the term coefficient
+            pauli_z = [Z(int(element[1:])) for element in pauli]
+            z_only_ham = SymbolicHamiltonian(term.coefficient * reduce(lambda x, y: x * y, pauli_z, 1.0))
+            # Can now apply expectation_from_samples directly
+            total += z_only_ham.expectation_from_samples(frequencies, qubit_map=qubits)
+        # Add the constant term if present. Note: Energies (in chemistry) are all real values
+        total += hamiltonian.constant.real
+        return total
+
+    # Expectation value from state vector simulation
+    result = circuit(nshots=1)
+    state_ket = result.state()
+    return hamiltonian.expectation(state_ket)

tests/test_expectation_samples.py

@@ -0,0 +1,69 @@
+"""
+Test expectation functionality
+"""
+
+# import numpy as np
+import pytest
+from qibo import Circuit, gates
+from qibo.hamiltonians import SymbolicHamiltonian
+from qibo.symbols import X, Z
+
+from qibochem.driver.molecule import Molecule
+from qibochem.measurement.expectation import expectation
+
+
+def test_expectation_z0():
+    """Test from_samples functionality of expectation function"""
+    hamiltonian = SymbolicHamiltonian(Z(0))
+    circuit = Circuit(2)
+    circuit.add(gates.X(0))
+    result = expectation(circuit, hamiltonian, from_samples=True)
+    assert pytest.approx(result) == -1.0
+
+
+def test_expectation_z0z1():
+    """Tests expectation_from_samples for diagonal Hamiltonians (only Z's)"""
+    hamiltonian = SymbolicHamiltonian(Z(0) * Z(1))
+    circuit = Circuit(2)
+    circuit.add(gates.X(0))
+    result = expectation(circuit, hamiltonian, from_samples=True)
+    assert pytest.approx(result) == -1.0
+
+
+def test_expectation_x0():
+    """Tests expectation_from_samples for Hamiltonians with X"""
+    hamiltonian = SymbolicHamiltonian(X(0))
+    circuit = Circuit(2)
+    circuit.add(gates.H(0))
+    result = expectation(circuit, hamiltonian, from_samples=True)
+    assert pytest.approx(result) == 1.0
+
+
+def test_expectation_x0_2():
+    """Test 2 of expectation_from_samples for Hamiltonians with X"""
+    hamiltonian = SymbolicHamiltonian(X(0))
+    circuit = Circuit(2)
+    result = expectation(circuit, hamiltonian, from_samples=True, n_shots=10000)
+    assert pytest.approx(result, abs=0.05) == 0.00
+
+
+def test_h2_hf_energy():
+    """Test HF energy of H2 molecule"""
+    # Hardcoded benchmark results
+    h2_ref_energy = -1.117349035
+
+    h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    try:
+        h2.run_pyscf()
+    except ModuleNotFoundError:
+        h2.run_psi4()
+
+    # JW-HF circuit
+    circuit = Circuit(4)
+    circuit.add(gates.X(_i) for _i in range(2))
+
+    # Molecular Hamiltonian and the HF expectation value
+    hamiltonian = h2.hamiltonian()
+    hf_energy = expectation(circuit, hamiltonian, from_samples=True, n_shots=10000)
+
+    assert h2_ref_energy == pytest.approx(hf_energy, abs=0.005)