diff --git a/doc/source/api-reference/ansatz.rst b/doc/source/api-reference/ansatz.rst index fb931f8..13f890a 100644 --- a/doc/source/api-reference/ansatz.rst +++ b/doc/source/api-reference/ansatz.rst @@ -2,12 +2,12 @@ Ansatz ====== -This section covers the API reference for various chemistry ansatzes. +This section covers the API reference for various chemistry circuit ansatzes. Hardware-efficient ------------------ -.. autofunction:: qibochem.ansatz.hardware_efficient.hea +.. autofunction:: qibochem.ansatz.hardware_efficient.he_circuit Hartree-Fock diff --git a/doc/source/getting-started/quickstart.rst b/doc/source/getting-started/quickstart.rst index a2ed28d..ee83b7b 100644 --- a/doc/source/getting-started/quickstart.rst +++ b/doc/source/getting-started/quickstart.rst @@ -15,8 +15,7 @@ Here is an example of building the UCCD ansatz with the H2 molecule to test your from qibo.models import VQE from qibochem.driver.molecule import Molecule - from qibochem.ansatz.hf_reference import hf_circuit - from qibochem.ansatz.ucc import ucc_circuit + from qibochem.ansatz import hf_circuit, ucc_circuit # Define the H2 molecule and obtain its 1-/2- electron integrals with PySCF h2 = Molecule([('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.7))]) diff --git a/doc/source/tutorials/ansatz.rst b/doc/source/tutorials/ansatz.rst index 3dc4d4f..b352f7c 100644 --- a/doc/source/tutorials/ansatz.rst +++ b/doc/source/tutorials/ansatz.rst @@ -16,17 +16,12 @@ For the H\ :sub:`2` case discussed in previous sections, a possible hardware eff .. code-block:: python - from qibo import Circuit + from qibochem.ansatz import he_circuit - from qibochem.ansatz import hardware_efficient - - nlayers = 1 nqubits = 4 - nfermions = 2 + nlayers = 1 - circuit = Circuit(4) - hardware_efficient_ansatz = hardware_efficient.hea(nlayers, nqubits) - circuit.add(hardware_efficient_ansatz) + circuit = he_circuit(nqubits, nlayers) print(circuit.draw()) .. code-block:: output @@ -43,11 +38,10 @@ The following example demonstrates how the energy of the H2 molecule is affected .. code-block:: python import numpy as np - from qibo import Circuit - from qibochem.driver.molecule import Molecule + from qibochem.driver import Molecule from qibochem.measurement.expectation import expectation - from qibochem.ansatz import hardware_efficient + from qibochem.ansatz import he_circuit mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.74804))]) mol.run_pyscf() @@ -57,10 +51,7 @@ The following example demonstrates how the energy of the H2 molecule is affected nlayers = 1 nqubits = mol.nso ntheta = 2 * nqubits * nlayers - hea_ansatz = hardware_efficient.hea(nlayers, nqubits) - - circuit = Circuit(nqubits) - circuit.add(hea_ansatz) + circuit = he_circuit(nqubits, nlayers) print("Energy expectation values for thetas: ") print("-----------------------------") @@ -116,8 +107,7 @@ An example of how to build a UCC doubles circuit ansatz for the :math:`H_2` mole .. code-block:: python from qibochem.driver.molecule import Molecule - from qibochem.ansatz.hf_reference import hf_circuit - from qibochem.ansatz.ucc import ucc_circuit + from qibochem.ansatz import hf_circuit, ucc_circuit mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.74804))]) mol.run_pyscf() diff --git a/examples/br_example.py b/examples/br_example.py index c0813b5..76bcde0 100644 --- a/examples/br_example.py +++ b/examples/br_example.py @@ -4,21 +4,15 @@ """ import numpy as np -from qibo import Circuit, gates, models -from qibo.optimizers import optimize -from scipy.optimize import minimize +from qibo.models import VQE -from qibochem.ansatz import basis_rotation -from qibochem.ansatz.hf_reference import hf_circuit -from qibochem.driver.molecule import Molecule +from qibochem.ansatz import basis_rotation, hf_circuit +from qibochem.driver import Molecule from qibochem.measurement.expectation import expectation # Define molecule and populate mol = Molecule(xyz_file="h3p.xyz") -try: - mol.run_pyscf() -except ModuleNotFoundError: - mol.run_psi4() +mol.run_pyscf() # Diagonalize H_core to use as the guess wave function @@ -59,16 +53,14 @@ def basis_rotation_circuit(mol, parameters=0.0): gate_layout = basis_rotation.basis_rotation_layout(nqubits) gate_list, ordered_angles = basis_rotation.basis_rotation_gates(gate_layout, gate_angles, theta) - circuit = Circuit(nqubits) - for _i in range(mol.nelec): - circuit.add(gates.X(_i)) + circuit = hf_circuit(nqubits, mol.nelec) circuit.add(gate_list) return circuit, gate_angles br_circuit, qubit_parameters = basis_rotation_circuit(mol, parameters=0.1) -vqe = models.VQE(br_circuit, hamiltonian) +vqe = VQE(br_circuit, hamiltonian) vqe_result = vqe.minimize(qubit_parameters) print(f" HF energy: {mol.e_hf:.8f}") diff --git a/examples/ucc_example1.py b/examples/ucc_example1.py index bb453f8..d8ffb9f 100644 --- a/examples/ucc_example1.py +++ b/examples/ucc_example1.py @@ -6,22 +6,18 @@ from qibo.optimizers import optimize from scipy.optimize import minimize -from qibochem.ansatz.hf_reference import hf_circuit -from qibochem.ansatz.ucc import ucc_circuit +from qibochem.ansatz import hf_circuit, ucc_circuit from qibochem.driver.molecule import Molecule from qibochem.measurement.expectation import expectation # Define molecule and populate mol = Molecule(xyz_file="lih.xyz") -try: - mol.run_pyscf() -except ModuleNotFoundError: - mol.run_psi4() +mol.run_pyscf() # Apply embedding and boson encoding mol.hf_embedding(active=[1, 2, 5]) -hamiltonian = mol.hamiltonian(oei=mol.embed_oei, tei=mol.embed_tei, constant=mol.inactive_energy) +hamiltonian = mol.hamiltonian() # Set parameters for the rest of the experiment n_qubits = mol.n_active_orbs @@ -56,11 +52,8 @@ n_unique_excitations = 5 # UCCSD: Circuit -all_coeffs = [] for _ex in excitations: - coeffs = [] - circuit += ucc_circuit(n_qubits, _ex, coeffs=coeffs) - all_coeffs.append(coeffs) + circuit += ucc_circuit(n_qubits, _ex) # Draw the circuit if interested print(circuit.draw()) @@ -71,19 +64,28 @@ def electronic_energy(parameters): r""" Loss function for the UCCSD ansatz """ - all_parameters = [] + coeff_dict = {1: (-1.0, 1.0), 2: (-0.25, 0.25, 0.25, 0.25, -0.25, -0.25, -0.25, 0.25)} - # UCC parameters - # Expand the parameters to match the total UCC ansatz manually - _ucc = parameters[:n_unique_excitations] + # Unique UCC parameters # Manually group the related excitations together - ucc_parameters = [_ucc[0], _ucc[1], _ucc[2], _ucc[2], _ucc[3], _ucc[3], _ucc[4], _ucc[4]] - # Need to iterate through each excitation this time - for _coeffs, _parameter in zip(all_coeffs, ucc_parameters): + ucc_parameters = [ + parameters[0], + parameters[1], + parameters[2], + parameters[2], + parameters[3], + parameters[3], + parameters[4], + parameters[4], + ] + all_parameters = [] + # Iterate through each excitation, with the list of coefficients dependent on whether S/D excitation + for _ex, _parameter in zip(excitations, ucc_parameters): + coeffs = coeff_dict[len(_ex) // 2] # Convert a single value to a array with dimension=n_param_gates - ucc_parameter = np.repeat(_parameter, len(_coeffs)) + ucc_parameter = np.repeat(_parameter, len(coeffs)) # Multiply by coeffs - ucc_parameter *= _coeffs + ucc_parameter *= coeffs all_parameters.append(ucc_parameter) # Flatten all_parameters into a single list to set the circuit parameters @@ -101,7 +103,7 @@ def electronic_energy(parameters): best, params, extra = optimize(electronic_energy, params) -print("\nResults using Qibo optimize:") +print("\nResults using Qibo optimize: (With HF embedding)") print(f"FCI energy: {fci_energy:.8f}") print(f" HF energy: {mol.e_hf:.8f} (Classical)") print(f"VQE energy: {best:.8f} (UCCSD ansatz)") @@ -115,7 +117,7 @@ def electronic_energy(parameters): result = minimize(electronic_energy, params) best, params = result.fun, result.x -print("\nResults using scipy.optimize:") +print("\nResults using scipy.optimize: (With HF embedding)") print(f"FCI energy: {fci_energy:.8f}") print(f" HF energy: {mol.e_hf:.8f} (Classical)") print(f"VQE energy: {best:.8f} (UCCSD ansatz)") @@ -123,7 +125,7 @@ def electronic_energy(parameters): # print("Optimized parameters:", params) -full_ham = mol.hamiltonian("f") +full_ham = mol.hamiltonian("f", oei=mol.oei, tei=mol.tei, constant=0.0) mol_fci_energy = mol.eigenvalues(full_ham)[0] print(f"\nFCI energy: {mol_fci_energy:.8f} (Full Hamiltonian)") diff --git a/examples/ucc_example2.py b/examples/ucc_example2.py deleted file mode 100644 index 4241876..0000000 --- a/examples/ucc_example2.py +++ /dev/null @@ -1,91 +0,0 @@ -""" -Sample script of how to build and run the UCCSD ansatz with HF embedding for LiH. -This example uses some utility functions to build the UCCSD ansatz; should be easier to use - -TODO: IN-PROGRESS!!! -""" - -import numpy as np -from qibo.optimizers import optimize -from scipy.optimize import minimize - -from qibochem.ansatz.hf_reference import hf_circuit -from qibochem.ansatz.ucc import ucc_ansatz -from qibochem.driver.molecule import Molecule -from qibochem.measurement.expectation import expectation - -# Define molecule and populate -mol = Molecule(xyz_file="lih.xyz") -try: - mol.run_pyscf() -except ModuleNotFoundError: - mol.run_psi4() - - -# Apply embedding and boson encoding -mol.hf_embedding(active=[1, 2, 5]) -hamiltonian = mol.hamiltonian(oei=mol.embed_oei, tei=mol.embed_tei, constant=mol.inactive_energy) - -# Set parameters for the rest of the experiment -n_qubits = mol.n_active_orbs -n_electrons = mol.n_active_e - -# Build circuit -circuit = hf_circuit(n_qubits, n_electrons) - -circuit += ucc_ansatz(mol) - -# Draw the circuit? -# print(circuit.draw()) -# print() - - -def electronic_energy(parameters): - """ - Loss function for the circuit ansatz - - TODO: IN-PROGRESS; Will probably be written up as another utility function in ucc.py? - """ - circuit.set_parameters(parameters) - - return expectation(circuit, hamiltonian) - - -# Reference energy -fci_energy = hamiltonian.eigenvalues()[0] - -n_parameters = len(circuit.get_parameters()) - -# Random initialization -params = np.random.rand(n_parameters) - -# NOTE: Circuit parameters not restricted to be equal within a single UCC excitation! - -best, params, extra = optimize(electronic_energy, params) - -print("\nResults using Qibo optimize:") -print(f"FCI energy: {fci_energy:.8f}") -print(f" HF energy: {mol.e_hf:.8f} (Classical)") -print(f"VQE energy: {best:.8f} (UCCSD ansatz)") -print() -print("Optimized parameters:", params) - - -# Scipy minimize -params = np.random.rand(n_parameters) - -result = minimize(electronic_energy, params) -best, params = result.fun, result.x - -print("\nResults using scipy.optimize:") -print(f"FCI energy: {fci_energy:.8f}") -print(f" HF energy: {mol.e_hf:.8f} (Classical)") -print(f"VQE energy: {best:.8f} (UCCSD ansatz)") -print() -print("Optimized parameters:", params) - - -full_ham = mol.hamiltonian("f") -mol_fci_energy = mol.eigenvalues(full_ham)[0] - -print(f"\nFCI energy: {mol_fci_energy:.8f} (Full Hamiltonian)") diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py index a94b14c..02cee9c 100644 --- a/src/qibochem/ansatz/__init__.py +++ b/src/qibochem/ansatz/__init__.py @@ -1,15 +1,12 @@ -from qibochem.ansatz.basis_rotation import ( - basis_rotation_gates, - basis_rotation_layout, - givens_qr_decompose, - unitary, -) -from qibochem.ansatz.hardware_efficient import hea -from qibochem.ansatz.hf_reference import bk_matrix, bk_matrix_power2, hf_circuit +from qibochem.ansatz.basis_rotation import basis_rotation_gates +from qibochem.ansatz.hardware_efficient import he_circuit +from qibochem.ansatz.hf_reference import hf_circuit from qibochem.ansatz.ucc import ( - expi_pauli, generate_excitations, mp2_amplitude, sort_excitations, + ucc_ansatz, ucc_circuit, ) + +# TODO: Probably can move some of the functions, e.g. generate_excitations/sort_excitations to a new 'util.py' diff --git a/src/qibochem/ansatz/hardware_efficient.py b/src/qibochem/ansatz/hardware_efficient.py index 89a5fef..306bd63 100644 --- a/src/qibochem/ansatz/hardware_efficient.py +++ b/src/qibochem/ansatz/hardware_efficient.py @@ -1,32 +1,38 @@ -from qibo import gates - - -def hea(n_layers, n_qubits, parameter_gates=["RY", "RZ"], coupling_gates="CZ"): - """ - 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 - n_qubits: Number of qubits in the quantum circuit - parameter_gates: List of single-qubit rotation gates to be used in the ansatz. The gates should be given as - strings representing valid ``Qibo`` one-qubit gates. Default: ``["RY", "RZ"]`` - coupling_gates: String representing the two-qubit entangling gate to be used in the ansatz; should be a - valid two-qubit ``Qibo`` gate. Default: ``"CZ"`` - - Returns: - ``list``: Gates corresponding to the hardware-efficient ansatz - """ - gate_list = [] - - for ilayer in range(n_layers): - for rgate in parameter_gates: - for iqubit in range(n_qubits): - gate_list.append(getattr(gates, rgate)(iqubit, theta=0)) - - # entanglement - for iqubit in range(n_qubits - 1): - gate_list.append(getattr(gates, coupling_gates)(iqubit, iqubit + 1)) - gate_list.append(getattr(gates, coupling_gates)(n_qubits - 1, 0)) - - return gate_list +""" +Hardware efficient ansatz, e.g. Kandala et al. (https://doi.org/10.1038/nature23879) +""" + +from qibo import Circuit, gates + + +def he_circuit(n_qubits, n_layers, parameter_gates=None, coupling_gates="CZ"): + """ + Builds the generalized hardware-efficient ansatz, in which the rotation and entangling gates used can be + chosen by the user + + Args: + n_qubits: Number of qubits in the quantum circuit + n_layers: Number of layers of rotation and entangling gates + parameter_gates: Iterable of single-qubit rotation gates to be used in the ansatz. The gates should be given as + strings representing valid ``Qibo`` one-qubit gates. Default: ``["RY", "RZ"]`` + coupling_gates: String representing the two-qubit entangling gate to be used in the ansatz; should be a + valid two-qubit ``Qibo`` gate. Default: ``"CZ"`` + + Returns: + Qibo ``Circuit``: Circuit corresponding to the hardware-efficient ansatz + """ + # Default variables + if parameter_gates is None: + parameter_gates = ["RY", "RZ"] + + circuit = Circuit(n_qubits) + + for _ in range(n_layers): + # Rotation gates + circuit.add(getattr(gates, rgate)(qubit, theta=0.0) for qubit in range(n_qubits) for rgate in parameter_gates) + + # Entanglement gates + circuit.add(getattr(gates, coupling_gates)(qubit, qubit + 1) for qubit in range(n_qubits - 1)) + circuit.add(getattr(gates, coupling_gates)(n_qubits - 1, 0)) + + return circuit diff --git a/src/qibochem/ansatz/ucc.py b/src/qibochem/ansatz/ucc.py index df38bbc..4fde184 100644 --- a/src/qibochem/ansatz/ucc.py +++ b/src/qibochem/ansatz/ucc.py @@ -4,62 +4,54 @@ import numpy as np import openfermion -from qibo import gates, models +from qibo import Circuit, gates from qibochem.ansatz.hf_reference import hf_circuit -def expi_pauli(n_qubits, theta, pauli_string): - r""" +def expi_pauli(n_qubits, pauli_string, theta): + """ Build circuit representing exp(i*theta*pauli_string) Args: n_qubits: No. of qubits in the quantum circuit - theta: parameter - pauli_string: OpenFermion QubitOperator object, e.g. :math: `X_0 Y_1 Z_2 X_4` + pauli_string: String in the format: ``"X0 Z1 Y3 X11"`` + theta: Real number Returns: circuit: Qibo Circuit object representing exp(i*theta*pauli_string) - coeff: Coefficient of theta. May be useful for VQE """ - # Unpack the dictionary from pauli_string.terms.items() into - # a (tuple of Pauli letters) and the coefficient of the pauli_string - [(p_letters, _coeff)] = pauli_string.terms.items() - - # _coeff is an imaginary number, i.e. exp(i\theta something) - # Convert to the real coefficient for multiplying theta - coeff = -2.0 * np.real(_coeff * -1.0j) - - # Generate the list of basis change gates using the p_letters list - basis_changes = [ - gates.H(_qubit) if _gate == "X" else gates.RX(_qubit, -0.5 * np.pi, trainable=False) - for _qubit, _gate in p_letters - if _gate != "Z" - ] + # Split pauli_string into the old p_letters format + pauli_ops = sorted(((int(_op[1:]), _op[0]) for _op in pauli_string.split()), key=lambda x: x[0]) + n_pauli_ops = len(pauli_ops) + + # Convert theta into a real number for applying with a RZ gate + rz_parameter = -2.0 * theta + + # Generate the list of basis change gates using the pauli_ops list + basis_changes = [] + for qubit, pauli_op in pauli_ops: + if pauli_op == "Y": + basis_changes.append(gates.S(qubit).dagger()) + if pauli_op not in ("I", "Z"): + basis_changes.append(gates.H(qubit)) # Build the circuit - circuit = models.Circuit(n_qubits) + circuit = Circuit(n_qubits) # 1. Change to X/Y where necessary - circuit.add(_gate for _gate in basis_changes) - # 2. Add CNOTs to all pairs of qubits in p_letters, starting from the last letter - circuit.add( - gates.CNOT(_qubit1, _qubit2) for (_qubit1, _g1), (_qubit2, _g2) in zip(p_letters[::-1], p_letters[::-1][1:]) - ) - # 3. Add RZ gate to last element of p_letters - circuit.add(gates.RZ(p_letters[0][0], coeff * theta)) - # 4. Add CNOTs to all pairs of qubits in p_letters - circuit.add(gates.CNOT(_qubit2, _qubit1) for (_qubit1, _g1), (_qubit2, _g2) in zip(p_letters, p_letters[1:])) + circuit.add(basis_changes) + # 2. Add CNOTs to all pairs of qubits in pauli_ops, starting from the last letter + circuit.add(gates.CNOT(pauli_ops[_i][0], pauli_ops[_i - 1][0]) for _i in range(n_pauli_ops - 1, 0, -1)) + # 3. Add RZ gate to last element of pauli_ops + circuit.add(gates.RZ(pauli_ops[0][0], rz_parameter)) + # 4. Add CNOTs to all pairs of qubits in pauli_ops + circuit.add(gates.CNOT(pauli_ops[_i + 1][0], pauli_ops[_i][0]) for _i in range(n_pauli_ops - 1)) # 3. Change back to the Z basis - # .dagger() doesn't keep trainable=False, so need to use a for loop - # circuit.add(_gate.dagger() for _gate in basis_changes) - for _gate in basis_changes: - gate = _gate.dagger() - gate.trainable = False - circuit.add(gate) - return circuit, coeff + circuit.add(_gate.dagger() for _gate in reversed(basis_changes)) + return circuit -def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map=None, coeffs=None): +def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map=None): r""" Circuit corresponding to the unitary coupled-cluster ansatz for a single excitation @@ -70,8 +62,6 @@ def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map theta: UCC parameter. Defaults to 0.0 trotter_steps: Number of Trotter steps; i.e. number of times the UCC ansatz is applied with ``theta`` = ``theta / trotter_steps``. Default: 1 ferm_qubit_map: Fermion-to-qubit transformation. Default is Jordan-Wigner (``jw``). - coeffs: List to hold the coefficients for the rotation parameter in each Pauli string. - May be useful in running the VQE. WARNING: Will be modified in this function Returns: Qibo ``Circuit``: Circuit corresponding to a single UCC excitation @@ -103,15 +93,18 @@ def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map # Apply the qubit_ucc_operator 'trotter_steps' times: assert trotter_steps > 0, f"{trotter_steps} must be > 0!" - circuit = models.Circuit(n_qubits) + circuit = Circuit(n_qubits) for _i in range(trotter_steps): # Use the get_operators() generator to get the list of excitation operators - for pauli_string in qubit_ucc_operator.get_operators(): - _circuit, coeff = expi_pauli(n_qubits, theta / trotter_steps, pauli_string) + for raw_pauli_string in qubit_ucc_operator.get_operators(): + # Convert each operator into a string and get the associated coefficient + ((pauli_ops, coeff),) = raw_pauli_string.terms.items() # Unpack the single-item dictionary + pauli_string = " ".join(f"{pauli_op[1]}{pauli_op[0]}" for pauli_op in pauli_ops) + # Build the circuit and add it on + _circuit = expi_pauli( + n_qubits, pauli_string, -1.0j * coeff * theta / trotter_steps + ) # Divide imag. coeff by 1.0j circuit += _circuit - if isinstance(coeffs, list): - coeffs.append(coeff) - return circuit @@ -175,7 +168,7 @@ def generate_excitations(order, excite_from, excite_to, conserve_spin=True): return [[]] # Generate all possible excitations first - from itertools import combinations + from itertools import combinations # pylint: disable=C0415 all_excitations = [ [*_from, *_to] for _from in combinations(excite_from, order) for _to in combinations(excite_to, order) @@ -219,8 +212,9 @@ def sort_excitations(excitations): prev = [] # No idea how I came up with this, but it seems to work for double excitations - # Default sorting is OK for single excitations - sorting_fn = lambda x: sum((order + 1 - _i) * abs(x[2 * _i + 1] // 2 - x[2 * _i] // 2) for _i in range(0, order)) + def sorting_fn(x): + # Default sorting is OK for single excitations + return sum((order + 1 - _i) * abs(x[2 * _i + 1] // 2 - x[2 * _i] // 2) for _i in range(0, order)) # Make a copy of the list of excitations, and use it populate a new list iteratively while copy_excitations: @@ -344,13 +338,7 @@ def ucc_ansatz( if include_hf: circuit = hf_circuit(n_orbs, n_elec, ferm_qubit_map=ferm_qubit_map) else: - circuit = models.Circuit(n_orbs) + circuit = Circuit(n_orbs) for excitation, theta in zip(excitations, thetas): - # coeffs = [] - circuit += ucc_circuit( - n_orbs, excitation, theta, trotter_steps=trotter_steps, ferm_qubit_map=ferm_qubit_map # , coeffs=coeffs) - ) - # if isinstance(all_coeffs, list): - # all_coeffs.append(np.array(coeffs)) - + circuit += ucc_circuit(n_orbs, excitation, theta, trotter_steps=trotter_steps, ferm_qubit_map=ferm_qubit_map) return circuit diff --git a/tests/test_hardware_efficient.py b/tests/test_hardware_efficient.py index 08ef08c..e9ad4c7 100644 --- a/tests/test_hardware_efficient.py +++ b/tests/test_hardware_efficient.py @@ -1,54 +1,57 @@ import numpy as np import pytest -from qibo import gates, models -from scipy.optimize import minimize +from qibo import Circuit, gates +from qibo.optimizers import optimize -from qibochem.ansatz import hardware_efficient -from qibochem.driver.molecule import Molecule +from qibochem.ansatz import he_circuit +from qibochem.driver import Molecule from qibochem.measurement.expectation import expectation -def test_hea_ansatz(): - mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.74804))]) - mol.run_pyscf() - mol_classical_hf_energy = mol.e_hf - mol_sym_ham = mol.hamiltonian("s") - - nlayers = 1 - nqubits = mol.nso - ntheta = 2 * nqubits * nlayers - theta = np.zeros(ntheta) - - hea_ansatz = hardware_efficient.hea(nlayers, nqubits) - qc = models.Circuit(nqubits) - qc.add(gates.X(_i) for _i in range(mol.nelec)) - qc.add(hea_ansatz) - qc.set_parameters(theta) - - hf_energy = expectation(qc, mol_sym_ham) - assert mol_classical_hf_energy == pytest.approx(hf_energy) - - -def test_vqe_hea_ansatz(): - def test_vqe_hea_ansatz_cost(parameters, circuit, hamiltonian): +def test_he_circuit(): + n_qubits = 4 + n_layers = 1 + rotation_gates = ["RX"] + entanglement_gate = "CNOT" + gate_list = [] + for _ in range(n_layers): + # Rotation gates + gate_list += [ + getattr(gates, rotation_gate)(_i, 0.0) for rotation_gate in rotation_gates for _i in range(n_qubits) + ] + # Entanglement gates + gate_list += [getattr(gates, entanglement_gate)(_i, _i + 1) for _i in range(n_qubits - 1)] + # Test function + test_circuit = he_circuit(n_qubits, n_layers, rotation_gates, entanglement_gate) + + # Check gates are correct + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(gate_list, list(test_circuit.queue)) + ) + # Check that only two parametrised gates + assert len(test_circuit.get_parameters()) == n_layers * n_qubits * len(rotation_gates) + + +def test_vqe_he_ansatz(): + # Loss function for VQE + def electronic_energy(parameters, circuit, hamiltonian): circuit.set_parameters(parameters) return expectation(circuit, hamiltonian) - mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.75))]) + mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) mol.run_pyscf() - mol_classical_hf_energy = mol.e_hf - mol_sym_ham = mol.hamiltonian("s") - - nlayers = 3 - nqubits = mol.nso - ntheta = 2 * nqubits * nlayers - theta = np.full(ntheta, np.pi / 4) - - hea_ansatz = hardware_efficient.hea(nlayers, nqubits, ["RY", "RZ"], "CNOT") - qc = models.Circuit(nqubits) - qc.add(gates.X(_i) for _i in range(mol.nelec)) - qc.add(hea_ansatz) - qc.set_parameters(theta) - - vqe_object = minimize(test_vqe_hea_ansatz_cost, theta, args=(qc, mol_sym_ham), method="Powell") - assert vqe_object.fun == pytest.approx(-1.1371170019838128) + hamiltonian = mol.hamiltonian() + + n_layers = 2 + n_qubits = mol.nso + rotation_gates = None + entanglement_gate = "CNOT" + circuit = he_circuit(n_qubits, n_layers, rotation_gates, entanglement_gate) + + n_parameters = len(circuit.get_parameters()) + thetas = np.full(n_parameters, 0.25 * np.pi) + best, params, extra = optimize(electronic_energy, thetas, args=(circuit, hamiltonian), method="BFGS") + + exact_result = -1.136189454 + assert best == pytest.approx(exact_result), f"{best} != exact_result (= {exact_result})" diff --git a/tests/test_ucc.py b/tests/test_ucc.py index b3d5581..56d70e1 100644 --- a/tests/test_ucc.py +++ b/tests/test_ucc.py @@ -2,14 +2,16 @@ Tests for the UCC ansatz and related functions """ -from functools import partial +from functools import reduce import numpy as np import pytest -from qibo import gates +from qibo import Circuit, gates, symbols +from qibo.hamiltonians import SymbolicHamiltonian from qibochem.ansatz.hf_reference import hf_circuit from qibochem.ansatz.ucc import ( + expi_pauli, generate_excitations, mp2_amplitude, sort_excitations, @@ -64,51 +66,83 @@ def test_mp2_amplitude_doubles(): @pytest.mark.parametrize( - "excitation,mapping,basis_rotations", + "pauli_string", [ - ([0, 2], None, ([("Y", 0), ("X", 2)], [("X", 0), ("Y", 2)])), # JW singles + "Z1", + "Z0 Z1", + "X0 X1", + "Y0 Y1", + ], +) +def test_expi_pauli(pauli_string): + n_qubits = 2 + theta = 0.1 + + # Build using exp(-i*theta*SymbolicHamiltonian) + pauli_ops = sorted(((int(_op[1]), _op[0]) for _op in pauli_string.split()), key=lambda x: x[0]) + control_circuit = Circuit(n_qubits) + pauli_term = SymbolicHamiltonian( + symbols.I(n_qubits - 1) + * reduce(lambda x, y: x * y, (getattr(symbols, pauli_op)(qubit) for qubit, pauli_op in pauli_ops)) + ) + control_circuit += pauli_term.circuit(-theta) + control_result = control_circuit(nshots=1) + control_state = control_result.state(True) + + test_circuit = expi_pauli(n_qubits, pauli_string, theta) + test_result = test_circuit(nshots=1) + test_state = test_result.state(True) + + assert np.allclose(control_state, test_state) + + +@pytest.mark.parametrize( + "excitation,mapping,basis_rotations,coeffs", + [ + ([0, 2], None, ("Y0 X2", "X0 Y2"), (0.5, -0.5)), # JW singles ( [0, 1, 2, 3], None, ( - [("X", 0), ("X", 1), ("Y", 2), ("X", 3)], - [("Y", 0), ("Y", 1), ("Y", 2), ("X", 3)], - [("Y", 0), ("X", 1), ("X", 2), ("X", 3)], - [("X", 0), ("Y", 1), ("X", 2), ("X", 3)], - [("Y", 0), ("X", 1), ("Y", 2), ("Y", 3)], - [("X", 0), ("Y", 1), ("Y", 2), ("Y", 3)], - [("X", 0), ("X", 1), ("X", 2), ("Y", 3)], - [("Y", 0), ("Y", 1), ("X", 2), ("Y", 3)], + "X0 X1 Y2 X3", + "Y0 Y1 Y2 X3", + "Y0 X1 X2 X3", + "X0 Y1 X2 X3", + "Y0 X1 Y2 Y3", + "X0 Y1 Y2 Y3", + "X0 X1 X2 Y3", + "Y0 Y1 X2 Y3", ), + (-0.25, 0.25, 0.25, 0.25, -0.25, -0.25, -0.25, 0.25), ), # JW doubles - ([0, 2], "bk", ([("X", 0), ("Y", 1), ("X", 2)], [("Y", 0), ("Y", 1), ("Y", 2)])), # BK singles + ([0, 2], "bk", ("X0 Y1 X2", "Y0 Y1 Y2"), (0.5, 0.5)), # BK singles ], ) -def test_ucc_circuit(excitation, mapping, basis_rotations): +def test_ucc_circuit(excitation, mapping, basis_rotations, coeffs): """Build a UCC circuit with only one excitation""" - gate_dict = {"X": gates.H, "Y": partial(gates.RX, theta=-0.5 * np.pi, trainable=False)} - # Build the list of basis rotation gates - basis_rotation_gates = [ - [gate_dict[_gate[0]](_gate[1]) for _gate in basis_rotation] for basis_rotation in basis_rotations - ] - # Build the CNOT cascade manually - cnot_cascade = [gates.CNOT(_i, _i - 1) for _i in range(excitation[-1], excitation[0], -1)] - cnot_cascade = cnot_cascade + [gates.RZ(excitation[0], 0.0)] - cnot_cascade = cnot_cascade + [gates.CNOT(_i + 1, _i) for _i in range(excitation[0], excitation[-1])] - - # Build control list of gates - nested_gate_list = [gate_list + cnot_cascade + gate_list for gate_list in basis_rotation_gates] - gate_list = [_gate for gate_list in nested_gate_list for _gate in gate_list] - - # Test ucc_function - circuit = ucc_circuit(4, excitation, ferm_qubit_map=mapping) - # Check gates are correct - assert all( - control.name == test.name and control.target_qubits == test.target_qubits - for control, test in zip(gate_list, list(circuit.queue)) - ) - # Check that only two parametrised gates - assert len(circuit.get_parameters()) == len(basis_rotations) + theta = 0.1 + n_qubits = 4 + + # Build the control array using SymbolicHamiltonian.circuit + # But need to multiply theta by some coefficient introduced by the fermion->qubit mapping + control_circuit = Circuit(n_qubits) + for coeff, basis_rotation in zip(coeffs, basis_rotations): + n_terms = len(basis_rotation) + pauli_term = SymbolicHamiltonian( + symbols.I(n_qubits - 1) + * reduce(lambda x, y: x * y, (getattr(symbols, _op)(int(qubit)) for _op, qubit in basis_rotation.split())) + ) + control_circuit += pauli_term.circuit(-coeff * theta) + control_result = control_circuit(nshots=1) + control_state = control_result.state(True) + # Test the ucc_circuit function + test_circuit = ucc_circuit(n_qubits, excitation, theta=theta, ferm_qubit_map=mapping) + test_result = test_circuit(nshots=1) + test_state = test_result.state(True) + assert np.allclose(control_state, test_state) + + # Check that number of parametrised gates matches + assert len(test_circuit.get_parameters()) == len(basis_rotations) def test_ucc_ferm_qubit_map_error(): @@ -117,23 +151,6 @@ def test_ucc_ferm_qubit_map_error(): ucc_circuit(2, [0, 1], ferm_qubit_map="Unknown") -def test_ucc_parameter_coefficients(): - """Coefficients used to multiply the parameters in the UCC circuit. Note: may change in future!""" - # UCC-JW singles - control_values = (-1.0, 1.0) - coeffs = [] - _circuit = ucc_circuit(2, [0, 1], coeffs=coeffs) - # Check that the signs of the coefficients have been saved - assert all(control == test for control, test in zip(control_values, coeffs)) - - # UCC-JW doubles - control_values = (-0.25, 0.25, 0.25, 0.25, -0.25, -0.25, -0.25, 0.25) - coeffs = [] - _circuit = ucc_circuit(4, [0, 1, 2, 3], coeffs=coeffs) - # Check that the signs of the coefficients have been saved - assert all(control == test for control, test in zip(control_values, coeffs)) - - def test_ucc_ansatz_h2(): """Test the default arguments of ucc_ansatz using H2""" mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])