diff --git a/doc/source/api-reference/ansatz.rst b/doc/source/api-reference/ansatz.rst index ece5596..79d71eb 100644 --- a/doc/source/api-reference/ansatz.rst +++ b/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 diff --git a/doc/source/api-reference/measurement.rst b/doc/source/api-reference/measurement.rst index 7123879..da2a1f6 100644 --- a/doc/source/api-reference/measurement.rst +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 7b40d51..16ccd25 100644 --- a/pyproject.toml +++ b/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", diff --git a/src/qibochem/ansatz/basis_rotation.py b/src/qibochem/ansatz/basis_rotation.py index 4a35dd7..82009ae 100644 --- a/src/qibochem/ansatz/basis_rotation.py +++ b/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 diff --git a/src/qibochem/ansatz/hardware_efficient.py b/src/qibochem/ansatz/hardware_efficient.py index c84c96a..5b13e17 100644 --- a/src/qibochem/ansatz/hardware_efficient.py +++ b/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 diff --git a/src/qibochem/measurement/basis_rotate.py b/src/qibochem/measurement/basis_rotate.py deleted file mode 100644 index 69ba37f..0000000 --- a/src/qibochem/measurement/basis_rotate.py +++ /dev/null @@ -1,27 +0,0 @@ -import qibo -from qibo import gates -from qibo.models import Circuit - - -def measure_rotate_basis(pauli_op, nqubits): - """ - rotate each qubit to computational basis (Z-basis) for measurement - input: - pauli_op: qubit operator - nqubits: number of qubits - output: - qibo circuit with rotations and measurement gates - """ - - mqc = Circuit(nqubits) - for p in pauli_op: - if p[1] == "Z": - mqc.add(gates.M(p[0])) - elif p[1] == "Y": - mqc.add(gates.S(p[0]).dagger()) - mqc.add(gates.H(p[0])) - mqc.add(gates.M(p[0])) - elif p[1] == "X": - mqc.add(gates.H(p[0])) - mqc.add(gates.M(p[0])) - return mqc diff --git a/src/qibochem/measurement/expectation.py b/src/qibochem/measurement/expectation.py index e0c5ff1..6d4d2b6 100644 --- a/src/qibochem/measurement/expectation.py +++ b/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) diff --git a/tests/test_expectation_samples.py b/tests/test_expectation_samples.py new file mode 100644 index 0000000..e863053 --- /dev/null +++ b/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)