From d63c22197a09a81e14b7455fbc1425848468672e Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Fri, 16 Feb 2024 03:15:34 +0000 Subject: [PATCH 01/32] Upload students' work as Python scripts --- examples/h2.xyz | 4 + examples/paper2.py | 58 +++++++++ examples/paper3_caleb.py | 122 +++++++++++++++++++ examples/paper3_conan.py | 109 +++++++++++++++++ examples/psr_paper2.py | 233 +++++++++++++++++++++++++++++++++++ examples/psr_paper3.py | 257 +++++++++++++++++++++++++++++++++++++++ 6 files changed, 783 insertions(+) create mode 100644 examples/h2.xyz create mode 100644 examples/paper2.py create mode 100644 examples/paper3_caleb.py create mode 100644 examples/paper3_conan.py create mode 100644 examples/psr_paper2.py create mode 100644 examples/psr_paper3.py diff --git a/examples/h2.xyz b/examples/h2.xyz new file mode 100644 index 0000000..6a4e902 --- /dev/null +++ b/examples/h2.xyz @@ -0,0 +1,4 @@ +2 + +H 0.0 0.0 0.0 +H 0.0 0.0 0.7 diff --git a/examples/paper2.py b/examples/paper2.py new file mode 100644 index 0000000..3159d6b --- /dev/null +++ b/examples/paper2.py @@ -0,0 +1,58 @@ +""" +A rough Python script attempting to implement the VQE circuit ansatz proposed in Gard et al. using Qibo. + +Reference paper: https://doi.org/10.1038/s41534-019-0240-1 + +Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Caleb Seow from +Eunoia Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for +Junior College students. +""" + +import numpy as np +from qibo import Circuit, gates + +from qibochem.driver.molecule import Molecule + +# Define molecule +mol = Molecule(xyz_file="h2.xyz") +mol.run_pyscf() + + +n = 4 # no of qubits +m = 2 # no of electrons +theta = 0.1 +phi = 0.1 + + +def special_R(a, theta, phi): + A = Circuit(n) + A.add(gates.RY(a, theta + np.pi / 2)) + A.add(gates.RZ(a, phi + np.pi)) + return A + + +print(special_R(1, 0.1, 0.1).invert().draw()) +print(special_R(1, 0.1, 0.1).draw()) + + +def A_gate(a): + A = Circuit(n) + A.add(gates.CNOT(a + 1, a)) + A += special_R(a + 1, theta, phi).invert() + A.add(gates.CNOT(a, a + 1)) + A += special_R(a + 1, theta, phi) + A.add(gates.CNOT(a + 1, a)) + return A + + +print(A_gate(0).draw()) + +Gate2 = Circuit(n) +Gate2.add(gates.X(1)) +Gate2.add(gates.X(2)) +for i in range(2): + Gate2 += A_gate(0) + Gate2 += A_gate(2) + Gate2 += A_gate(1) + +print(Gate2.draw()) diff --git a/examples/paper3_caleb.py b/examples/paper3_caleb.py new file mode 100644 index 0000000..122bc37 --- /dev/null +++ b/examples/paper3_caleb.py @@ -0,0 +1,122 @@ +""" +A rough Python script attempting to implement the VQE circuit ansatz proposed in Arrazola et al. using Qibo. + +Reference paper: https://doi.org/10.22331/q-2022-06-20-742 + +Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Caleb Seow from +Eunoia Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for +Junior College students. +""" + +import numpy as np +from qibo import Circuit, gates, models + +from qibochem.ansatz.hf_reference import hf_circuit +from qibochem.driver.molecule import Molecule +from qibochem.measurement.expectation import expectation + +mol = Molecule(xyz_file="h2.xyz") +mol.run_pyscf() +n = 4 + +# Single Excitation + +g1 = hf_circuit(mol.nso, 2) + +theta = 0.1 + + +def single_ex(a, b): + g = Circuit(4) + g.add(gates.CNOT(a, b)) + g.add(gates.RY(a, theta)) + g.add(gates.CNOT(b, a)) + g.add(gates.RY(a, -theta)) + g.add(gates.CNOT(b, a)) + g.add(gates.CNOT(a, b)) + return g + + +for i in range(2): + for j in range(2, 4): + if (i + j) % 2 == 0: + g1 += single_ex(i, j) + +print(g1.summary()) +print(g1.draw()) + +theta = 0.1 + + +def double_ex(): + g = Circuit(4) + g.add(gates.CNOT(2, 3)) + g.add(gates.CNOT(0, 2)) + g.add(gates.H(0)) + g.add(gates.H(3)) + g.add(gates.CNOT(0, 1)) # 2 simultaneous CNOTs? + g.add(gates.CNOT(2, 3)) + g.add(gates.RY(0, -theta)) + g.add(gates.RY(1, theta)) + g.add(gates.CNOT(0, 3)) + g.add(gates.H(3)) + g.add(gates.CNOT(3, 1)) + g.add(gates.RY(0, -theta)) + g.add(gates.RY(1, theta)) + g.add(gates.CNOT(2, 1)) + g.add(gates.CNOT(2, 0)) + g.add(gates.RY(0, theta)) + g.add(gates.RY(1, -theta)) + g.add(gates.CNOT(3, 1)) + g.add(gates.H(3)) + g.add(gates.CNOT(0, 3)) + g.add(gates.RY(0, theta)) + g.add(gates.RY(1, -theta)) + g.add(gates.CNOT(0, 1)) + g.add(gates.CNOT(2, 0)) + g.add(gates.H(0)) + g.add(gates.H(3)) + g.add(gates.CNOT(0, 2)) + g.add(gates.CNOT(2, 3)) + return g + + +g1 += double_ex() + +print(g1.draw()) + +hamiltonian = mol.hamiltonian() +# parameter = float(input("theta = ")) +parameter = 0.1 +params = np.zeros(len(g1.get_parameters())) + parameter +g1.set_parameters(params) + +energy = expectation(g1, hamiltonian) +print(f"Expectation: {energy}") + +vqe = models.VQE(g1, hamiltonian) + +best, params, extra = vqe.minimize(params, method="BFGS", compile=False) +# Methods: BFGS, COBYLA, Powell + +print(f"Energy: {best}") +exact_result = hamiltonian.eigenvalues()[0] +print(f"Exact result: {exact_result}") + +hamiltonian = mol.hamiltonian() +# parameter = float(input("theta = ")) +parameter = 0.1 +params = np.zeros(len(g1.get_parameters())) + parameter +g1.set_parameters(params) + +energy = expectation(g1, hamiltonian) +print(f"Expectation: {energy}") + +vqe = models.VQE(g1, hamiltonian) + +best, params, extra = vqe.minimize(params, method="BFGS", compile=False) +# Methods: BFGS, COBYLA, Powell + +print(f"Energy: {best}") + +print(f"Exact result: {exact_result}") diff --git a/examples/paper3_conan.py b/examples/paper3_conan.py new file mode 100644 index 0000000..4ebcd23 --- /dev/null +++ b/examples/paper3_conan.py @@ -0,0 +1,109 @@ +""" +A rough Python script attempting to implement the VQE circuit ansatz proposed in Arrazola et al. using Qibo. + +Reference paper: https://doi.org/10.22331/q-2022-06-20-742 + +Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Conan Tan from +National Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for +Junior College students. +""" + +import numpy as np +from qibo import Circuit, gates, models + +from qibochem.ansatz.hf_reference import hf_circuit +from qibochem.driver.molecule import Molecule +from qibochem.measurement.expectation import expectation + +# Excitation function + + +def checkexcitation(elec, orb): + s_excitations = [(i, j) for i in range(elec) for j in range(elec, orb) if (i + j) % 2 == 0] + print(s_excitations) + + d_excitations = [ + (0, 1, k, l) + for k in range(elec, orb) + for l in range(k + 1, orb) + if (1 + k + l) % 2 == 0 and ((k % 2) + (l % 2)) == 1 + ] + print(d_excitations) + return s_excitations, d_excitations + + +# H2 + +# Define molecule and populate +mol = Molecule(xyz_file="h2.xyz") +mol.run_pyscf() +n_qubits = mol.nso +n_electrons = 2 +hamiltonian = mol.hamiltonian() +s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) + + +## Circuit construction +c = hf_circuit(n_qubits, n_electrons) + + +def build(c, singleex, doubleex, x): + for qa, qb in singleex: + sc = Circuit(n_qubits) + sc.add(gates.CNOT(qa, qb)) + sc.add(gates.RY(qa, theta=x / 2)) + sc.add(gates.CNOT(qb, qa)) + sc.add(gates.RY(qa, theta=-x / 2)) + sc.add(gates.CNOT(qb, qa)) + sc.add(gates.CNOT(qa, qb)) + c += sc + # for _i, parameter in enumerate(c.get_parameters()): + # gradient = parameter_shift(c, hamiltonian, parameter_index=_i) + # print(f"Excitation {qa, qb} => Gradient: {gradient}") + # if np.abs(gradient) > 1e-10: + # gradients[(qa, qb)] = gradient + # break + for qa, qb, qc, qd in doubleex: + dc = Circuit(n_qubits) + dc.add(gates.CNOT(qc, qd)) + dc.add(gates.CNOT(qa, qc)) + dc.add(gates.H(qa)) + dc.add(gates.H(qd)) + dc.add(gates.CNOT(qa, qb)) + dc.add(gates.CNOT(qc, qd)) + dc.add(gates.RY(qa, theta=-x / 8)) + dc.add(gates.RY(qb, theta=x / 8)) + dc.add(gates.CNOT(qa, qd)) + dc.add(gates.H(qd)) + dc.add(gates.CNOT(qd, qb)) + dc.add(gates.RY(qa, theta=-x / 8)) + dc.add(gates.RY(qb, theta=x / 8)) + dc.add(gates.CNOT(qc, qb)) + dc.add(gates.CNOT(qc, qa)) + dc.add(gates.RY(qa, theta=x / 8)) + dc.add(gates.RY(qb, theta=-x / 8)) + dc.add(gates.CNOT(qd, qb)) + dc.add(gates.H(qd)) + dc.add(gates.CNOT(qa, qd)) + dc.add(gates.RY(qa, theta=x / 8)) + dc.add(gates.RY(qb, theta=-x / 8)) + dc.add(gates.CNOT(qa, qb)) + dc.add(gates.CNOT(qc, qa)) + dc.add(gates.H(qa)) + dc.add(gates.H(qd)) + dc.add(gates.CNOT(qa, qc)) + dc.add(gates.CNOT(qc, qd)) + c += dc + # for _i, parameter in enumerate(c.get_parameters()): + # gradient = parameter_shift(c, hamiltonian, parameter_index=_i) + # print(f"Excitation {qa, qb, qc, qd} => Gradient: {gradient}") + # if np.abs(gradient) > 1e-10: + # gradients[(qa, qb, qc, qd)] = gradient + # break + return c + + +print(list(enumerate(c.get_parameters()))) +c = build(c, s_excitations, d_excitations, 0.1) +print(c.draw()) +print(c.summary()) diff --git a/examples/psr_paper2.py b/examples/psr_paper2.py new file mode 100644 index 0000000..856856e --- /dev/null +++ b/examples/psr_paper2.py @@ -0,0 +1,233 @@ +""" +This script is a continuation of 'paper2.py'. After implementing the circuit ansatz, this script attempts to examine the +use of the Parameter Shift Rule implemented in Qibo for chemistry applications. In addition, simulations of LiH are used +for some benchmarking. + +Reference paper: https://doi.org/10.1038/s41534-019-0240-1 + +Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Caleb Seow from +Eunoia Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for +Junior College students. +""" + +import math + +import matplotlib.pyplot as plt +import numpy as np +from qibo import Circuit, gates, models +from qibo.derivative import parameter_shift +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.driver.molecule import Molecule +from qibochem.measurement.expectation import expectation + +# mol = Molecule(xyz_file="h2.xyz") +h2 = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) +h2.run_pyscf() + +h2.hf_embedding(active=[1, 2, 5]) + +print(f"# electrons: {h2.n_active_e}") +print(f"# qubits: {h2.n_active_orbs}") + +# circuit = hf_circuit(h2.n_active_orbs, h2.n_active_e) +# print(circuit.draw()) +# print(h2.eps) +# quit() + +hamiltonian = h2.hamiltonian(oei=h2.embed_oei, tei=h2.embed_tei, constant=h2.inactive_energy) + +# Exact_result = hamiltonian.eigenvalues()[0] +# print(f"Exact result: {Exact_result}") +# quit() + +theta, phi = 0.0, 0.0 + + +def special_R(a, theta, phi): + A = Circuit(h2.n_active_orbs) + A.add(gates.RY(a, theta)) + A.add(gates.RZ(a, phi)) + return A + + +def A_gate(a): + A = Circuit(h2.n_active_orbs) + A.add(gates.CNOT(a + 1, a)) + A += special_R(a + 1, theta, phi).invert() + A.add(gates.CNOT(a, a + 1)) + A += special_R(a + 1, theta, phi) + A.add(gates.CNOT(a + 1, a)) + return A + + +# print(circuit.draw()) +# print(A_gate(0).get_parameters()) # 4 gradients + +# initial state +if h2.n_active_e <= math.ceil(h2.n_active_orbs): + Xpositions = [2 * i for i in range(h2.n_active_e)] +else: + Xpositions = [2 * i for i in range(math.ceil(h2.n_active_orbs / 2))] + Xpositions += [2 * j + 1 for j in range(h2.n_active_e - math.ceil(h2.n_active_orbs / 2))] + Xpositions.sort() + + +# setting the circuit up +def setup_electrons(g2): + for i in Xpositions: + g2.add(gates.X(i)) + + +Xpositions = [0, 2] + +print("Xpositions: ", Xpositions) + +# Xpositions is like a starting point +no_of_gates = 0 +var_list = [] +master_list = [] +ref_list = [] +popping_list = [] + +# initialising +for i in range(h2.n_active_orbs): + if i + 1 not in Xpositions: + if i != h2.n_active_orbs - 1: + ref_list.append((i, i + 1)) + master_list.append((i, i + 1)) + no_of_gates += 1 + +print(f"initial master list: {master_list}") + + +# Adding 1 layer based on the prev layer +while no_of_gates < math.comb(h2.n_active_orbs, h2.n_active_e): + + for i, j in ref_list: + var_list.append((i - 1, i)) + var_list.append((j, j + 1)) + + var_list.append((0, 0)) # dummy tuple that won't appear naturally + + for i in range(len(var_list) - 1): # 3 conditions to pop + if var_list[i][0] < 0: + popping_list.append(i) + elif var_list[i][1] > h2.n_active_orbs - 1: + popping_list.append(i) + elif var_list[i] == var_list[i + 1]: + popping_list.append(i + 1) + + popping_list.reverse() # so that index won't change after every pop + + for i in popping_list: + var_list.pop(i) + + var_list.pop() # removing the (0,0) tuple + + for i in var_list: + master_list.append(i) + no_of_gates += 1 + + ref_list.clear() + for i in var_list: + ref_list.append(i) + + popping_list.clear() + var_list.clear() + +print(f"final master list: {master_list}") + +A_gate_indexes = master_list[: math.comb(h2.n_active_orbs, h2.n_active_e)] # shouldnt be needed but jic?? + +print(f"A gate indexes: {A_gate_indexes}") + +g2 = Circuit(h2.n_active_orbs) + +setup_electrons(g2) + +for i, j in A_gate_indexes: + g2 += A_gate(i) + +print(g2.draw()) + +t, p = 0.1, 0.1 # theta, phi +params = [] + +# qn: for the dagger do we want the changed params too? +# A_gate_param = [-p-np.pi,-t-np.pi/2,t+np.pi/2,p+np.pi] +# A_gate_param = [p+np.pi,t+np.pi/2,t+np.pi/2,p+np.pi] +A_gate_param = [-t, -p, p, t] + +for i in range(len(A_gate_indexes)): + params += A_gate_param + +g2.set_parameters(params) + +# gradlist = [] +# +# gradients = {} +# +# print(len(g2.get_parameters())) +# +# thetalist = [f"theta {i}" for i in range(len(g2.get_parameters()))] +# +# for _i, parameter in enumerate(g2.get_parameters()): +# gradient = parameter_shift(g2, hamiltonian, parameter_index=_i) +# gradlist.append(gradient) +# +# for i,j in enumerate(thetalist): +# if np.abs(gradlist[i]) > 1e-10: +# gradients[j] = gradient + +# print(gradients) + +energy = expectation(g2, hamiltonian) +print(f"Expectation: {energy}") + +# vqe = models.VQE(g2, hamiltonian) + +# Methods: BFGS, COBYLA, Powell +# best_BFGS, params_BFGS, extra_BFGS = vqe.minimize(params, method='BFGS', compile=False) +# print(f"Energy BFGS: {best_BFGS}") +# +# best_COBYLA, params_COBYLA, extra_COBYLA = vqe.minimize(params, method='COBYLA', compile=False) +# print(f"Energy COBYLA: {best_COBYLA}") +# +# best_Powell, params_Powell, extra_Powell = vqe.minimize(params, method='Powell', compile=False) +# print(f"Energy Powell: {best_Powell}") + +# Exact_result = hamiltonian.eigenvalues()[0] +# print(f"Exact result: {Exact_result}") + + +def no_restriction_electronic_energy(parameters): + g2.set_parameters(parameters) + return expectation(g2, hamiltonian) + + +# Combine the gradients of each parameterised gate into a single array (vector): +def gradient_no_restriction(_thetas): + g2.set_parameters(_thetas) + return np.array([parameter_shift(g2, hamiltonian, parameter_index=_i) for _i in range(len(_thetas))]) + + +# Start from zero +thetas = np.zeros(len(g2.get_parameters())) + +best, params, extra = optimize(no_restriction_electronic_energy, thetas, method="SLSQP", jac=gradient_no_restriction) +# best, params, extra = vqe.minimize(thetas, method='BFGS', jac=ucc_gradient_no_restriction) + +print(f"VQE result: {best}") +print(f"Optimized parameter: {params}") +print(f"Number of steps: {extra.nfev}") +print() + +# COBYLA +best, params, extra = optimize(no_restriction_electronic_energy, thetas, method="SLSQP") +print(f"VQE result: {best}") +print(f"Optimized parameter: {params}") +print(f"Number of steps: {extra.nfev}") diff --git a/examples/psr_paper3.py b/examples/psr_paper3.py new file mode 100644 index 0000000..4ae8f8c --- /dev/null +++ b/examples/psr_paper3.py @@ -0,0 +1,257 @@ +""" +This script is a continuation of 'paper3_conan.py'. After implementing the circuit ansatz, this script attempts to +examine the use of the Parameter Shift Rule implemented in Qibo for chemistry applications. In addition, simulations of +LiH are used for some benchmarking. + +Reference paper: https://doi.org/10.22331/q-2022-06-20-742 + +Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Conan Tan from +National Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for +Junior College students. +""" + +import numpy as np +from qibo import Circuit, gates +from qibo.derivative import parameter_shift +from qibo.optimizers import optimize + +from qibochem.ansatz.hf_reference import hf_circuit +from qibochem.driver.molecule import Molecule +from qibochem.measurement.expectation import expectation + + +# Excitation function +def checkexcitation(elec, orb): + s_excitations = [(i, j) for i in range(elec) for j in range(elec, orb) if (i + j) % 2 == 0] + print(s_excitations) + + d_excitations = [ + (0, 1, k, l) + for k in range(elec, orb) + for l in range(k + 1, orb) + if (1 + k + l) % 2 == 0 and ((k % 2) + (l % 2)) == 1 + ] + print(d_excitations) + return s_excitations, d_excitations + + +# H2 +# Define molecule and populate +mol = Molecule(xyz_file="h2.xyz") +mol.run_pyscf() +n_qubits = mol.nso +n_electrons = mol.nelec +hamiltonian = mol.hamiltonian() +s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) +gradients = {} + +## Circuit construction +c = hf_circuit(n_qubits, n_electrons) + + +def build(c, singleex, doubleex, x): + for qa, qb in singleex: + sc = Circuit(n_qubits) + sc.add(gates.CNOT(qa, qb)) + sc.add(gates.RY(qa, theta=x / 2)) + sc.add(gates.CNOT(qb, qa)) + sc.add(gates.RY(qa, theta=-x / 2)) + sc.add(gates.CNOT(qb, qa)) + sc.add(gates.CNOT(qa, qb)) + c += sc + for _i, parameter in enumerate(c.get_parameters()): + gradient = parameter_shift(c, hamiltonian, parameter_index=_i) + print(f"Excitation {qa, qb} => Gradient: {gradient}") + if np.abs(gradient) > 1e-10: + gradients[(qa, qb)] = gradient + # break + + # for qa, qb, qc, qd in doubleex: + # dc = Circuit(n_qubits) + # dc.add(gates.CNOT(qc, qd)) + # dc.add(gates.CNOT(qa, qc)) + # dc.add(gates.H(qa)) + # dc.add(gates.H(qd)) + # dc.add(gates.CNOT(qa, qb)) + # dc.add(gates.CNOT(qc, qd)) + # dc.add(gates.RY(qa, theta=-x/8)) + # dc.add(gates.RY(qb, theta=x/8)) + # dc.add(gates.CNOT(qa, qd)) + # dc.add(gates.H(qd)) + # dc.add(gates.CNOT(qd, qb)) + # dc.add(gates.RY(qa, theta=-x/8)) + # dc.add(gates.RY(qb, theta=x/8)) + # dc.add(gates.CNOT(qc, qb)) + # dc.add(gates.CNOT(qc, qa)) + # dc.add(gates.RY(qa, theta=x/8)) + # dc.add(gates.RY(qb, theta=-x/8)) + # dc.add(gates.CNOT(qd, qb)) + # dc.add(gates.H(qd)) + # dc.add(gates.CNOT(qa, qd)) + # dc.add(gates.RY(qa, theta=x/8)) + # dc.add(gates.RY(qb, theta=-x/8)) + # dc.add(gates.CNOT(qa, qb)) + # dc.add(gates.CNOT(qc, qa)) + # dc.add(gates.H(qa)) + # dc.add(gates.H(qd)) + # dc.add(gates.CNOT(qa, qc)) + # dc.add(gates.CNOT(qc, qd)) + # c += dc + # for _i, parameter in enumerate(c.get_parameters()): + # gradient = parameter_shift(c, hamiltonian, parameter_index=_i) + # print(f"Excitation {qa, qb, qc, qd} => Gradient: {gradient}") + # if np.abs(gradient) > 1e-10: + # gradients[(qa, qb, qc, qd)] = gradient + # break + return c + + +# checkexcitation(sum(n_electrons), n_qubits) +print(list(enumerate(c.get_parameters()))) +c = build(c, s_excitations, d_excitations, 0.0) +print(c.draw()) +print(c.summary()) + +for excitation in sorted(gradients, key=lambda x: np.abs(gradients[x]), reverse=True): + print(f"Excitation {excitation} => Gradient: {gradients[excitation]}") + +params = [] +s_coeffs = [1 / 2, -1 / 2] +d_coeffs = [-1 / 8, 1 / 8, -1 / 8, 1 / 8, 1 / 8, -1 / 8, 1 / 8, -1 / 8] + +for i in s_excitations: + params += s_coeffs + +for i in d_excitations: + params += d_coeffs + +## Restricted +x = [] + + +def electronic_energy(x): + x *= params + c.set_parameters(x) + return expectation(c, hamiltonian) + + +y = [] + + +def gradient_restriction(y): + y *= params + c.set_parameters(y) + return np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))]) + + +print(np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))])) + + +## Non-Restricted +def nr_electronic_energy(parameter): + c.set_parameters(parameter) + return expectation(c, hamiltonian) + + +y = [] + + +def gradient_no_restriction(y): + return np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))]) + + +print(np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))])) + + +## Energy Calculation +methods = ["SLSQP", "BFGS", "COBYLA", "Powell"] + +for m in methods: + theta = np.zeros(len(c.get_parameters())) + if m in ["BFGS", "SLSQP"]: + # Gradients for BFGS and SLSQP + best, params, extra = optimize(electronic_energy, theta, method=m, jac=gradient_restriction) + else: + best, params, extra = optimize(electronic_energy, theta, method=m) + + print(f"Energy {m}: {best}") + print(f"Optimized parameter: {params}") + print(f"Number of steps: {extra.nfev}") + +print(f"Exact result: {hamiltonian.eigenvalues()[0]}") + + +# LiH +mol = Molecule(xyz_file="lih.xyz") +mol.run_pyscf() +mol.hf_embedding(active=[1, 2, 5]) +hamiltonian = mol.hamiltonian(oei=mol.embed_oei, tei=mol.embed_tei, constant=mol.inactive_energy) +n_qubits = mol.n_active_orbs +n_electrons = mol.n_active_e +s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) +gradients = {} + +## Circuit construction + +c = hf_circuit(n_qubits, n_electrons) +s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) +c = build(c, s_excitations, d_excitations, 0.1) +print(c.draw()) +print(c.summary()) + +print(f"\nInitial number of excitation: {len(s_excitations + d_excitations)}") +print(f"Final number of excitations: {len(gradients)}") + +for excitation in sorted(gradients, key=lambda x: np.abs(gradients[x]), reverse=True): + print(f"Excitation {excitation} => Gradient: {gradients[excitation]}") + +## Calculation of Hamiltonian + +param_range = np.arange(*map(float, input("Enter space-seperated parameter range, and step increment = ").split())) +for x in param_range: + print(f"Theta: {round(x, 5)} => Energy expectation value: {electronic_energy(x)}") + +## Restricted +x = [] + + +def electronic_energy(x): + x *= params + c.set_parameters(x) + return expectation(c, hamiltonian) + + +# param_range = np.arange(*map(float, input("Enter space-seperated parameter range, and step increment = ").split())) +# for x in param_range: +# print(f"Theta: {round(x, 5)} => Energy expectation value: {electronic_energy(x)}") + +# print(params) + +y = [] + + +def gradient_restriction(y): + y *= params + c.set_parameters(y) + return np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))]) + + +print(np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))])) + +## Energy Calculation + +methods = ["SLSQP", "BFGS", "COBYLA", "Powell"] + +for m in methods: + theta = np.zeros(len(c.get_parameters())) + if m in ["BFGS", "SLSQP"]: + # Gradients for BFGS and SLSQP + best, params, extra = optimize(electronic_energy, theta, method=m, jac=gradient_restriction) + else: + best, params, extra = optimize(electronic_energy, theta, method=m) + + print(f"Energy {m}: {best}") + print(f"Optimized parameter: {params}") + print(f"Number of steps: {extra.nfev}") + +print(f"Exact result: {hamiltonian.eigenvalues()[0]}") From 6b9862e8b31f7143f5baaba6d5462e802e301090 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 16 May 2024 09:48:58 +0000 Subject: [PATCH 02/32] Initial outline of paper2 circuit ansatz --- examples/psr_paper2.py | 10 ++++--- src/qibochem/ansatz/symmetry.py | 52 +++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 4 deletions(-) create mode 100644 src/qibochem/ansatz/symmetry.py diff --git a/examples/psr_paper2.py b/examples/psr_paper2.py index 856856e..37ddbb4 100644 --- a/examples/psr_paper2.py +++ b/examples/psr_paper2.py @@ -19,10 +19,9 @@ 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.driver.molecule import Molecule -from qibochem.measurement.expectation import expectation +from qibochem.ansatz import hf_circuit, ucc_circuit +from qibochem.driver import Molecule +from qibochem.measurement import expectation # mol = Molecule(xyz_file="h2.xyz") h2 = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) @@ -75,6 +74,9 @@ def A_gate(a): Xpositions += [2 * j + 1 for j in range(h2.n_active_e - math.ceil(h2.n_active_orbs / 2))] Xpositions.sort() +print(Xpositions) +quit() + # setting the circuit up def setup_electrons(g2): diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py new file mode 100644 index 0000000..2ab48d8 --- /dev/null +++ b/src/qibochem/ansatz/symmetry.py @@ -0,0 +1,52 @@ +""" +Symmetry-Preserving circuit ansatz from Gard et al. Reference: https://doi.org/10.1038/s41534-019-0240-1 +""" + +import numpy as np +from qibo import Circuit, gates + +# Helper functions + + +def a_gate(qubit1, qubit2): + """ + Returns the list of elementary gates corresponding to the 'A' gate as defined in the paper, acting on qubit1 and qubit2 + """ + result = [] + result.append(gates.CNOT(qubit2, qubit1)) + result += [gates.RY(qubit2, 0.0), gates.RZ(qubit2, 0.0)] + result.append(gates.CNOT(qubit1, qubit2)) + result += [gates.RZ(qubit2, 0.0), gates.RY(qubit2, 0.0)] + result.append(gates.CNOT(qubit2, qubit1)) + return result + + +# Main function +def symm_preserving_circuit(n_qubits, n_electrons): + """ + Symmetry-preserving circuit ansatz + + Args: + n_qubits: Number of qubits in the quantum circuit + n_electrons: Number of electrons in the molecular system + + Returns: + Qibo ``Circuit``: Circuit ansatz + """ + circuit = Circuit(n_qubits) + circuit.add(gates.X(2 * _i) for _i in range(n_electrons)) + + a_gate_qubits = [] # Generate the list of qubits pairs for adding A gates + + a_gate_qubits = [(0, 1)] + + a_gates = [a_gate(qubit1, qubit2) for qubit1, qubit2 in a_gate_qubits] + circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate) # Unpack the nested list + + return circuit + + +n_qubits = 4 +n_electrons = 2 +circuit = symm_preserving_circuit(n_qubits, n_electrons) +print(circuit.draw()) From b0d94e4fe9dbccdabb36f762bf4d7842e0e84ffe Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 27 Jun 2024 09:14:31 +0000 Subject: [PATCH 03/32] Finish first draft of symmetry-preserving ansatz Now testing --- src/qibochem/ansatz/symmetry.py | 88 ++++++++++++++++++++++++--------- 1 file changed, 65 insertions(+), 23 deletions(-) diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py index 2ab48d8..a753721 100644 --- a/src/qibochem/ansatz/symmetry.py +++ b/src/qibochem/ansatz/symmetry.py @@ -8,17 +8,68 @@ # Helper functions -def a_gate(qubit1, qubit2): +def a_gate(qubit1, qubit2, theta=None, phi=None): """ - Returns the list of elementary gates corresponding to the 'A' gate as defined in the paper, acting on qubit1 and qubit2 + Decomposition of the 'A' gate as defined in the paper, acting on qubit1 and qubit2. 'A' corresponds to the following + unitary matrix: + + A(\\theta, \\phi) = + \\begin{pmatrix} + 1 & 0 & 0 & 0 \\\\ + 0 & \\cos \\theta & e^{i \\phi} \\sin \\theta & 0 \\\\ + 0 & e^{-i \\phi} \\sin \\theta & -\\cos \\theta & 0 \\\\ + 0 & 0 & 0 & 1 + \\end{pmatrix} + + Args: + qubit1 (int): Index of the first qubit + qubit2 (int): Index of the second qubit + theta (float): First rotation angle. Default: 0.0 + phi (float): Second rotation angle. Default: 0.0 + + Returns: + (list): List of gates representing the decomposition of the 'A' gate + """ + if theta is None: + theta = 0.0 + if phi is None: + phi = 0.0 + + # R(theta, phi) = R_z (phi + pi) R_y (theta + 0.5*pi) + r_and_cnot = [gates.RZ(qubit2, phi + np.pi), gates.RY(qubit2, theta + 0.5 * np.pi), gates.CNOT(qubit2, qubit1)] + return r_and_cnot[::-1] + [gates.CNOT(qubit1, qubit2)] + r_and_cnot + + +def x_gate_indices(n_qubits, n_electrons): + """Obtain the qubit indices for X gates to be added to the circuit""" + indices = [2 * _i for _i in range(0, min(n_electrons, n_qubits // 2))] + if n_electrons > n_qubits // 2: + indices += [2 * _i + 1 for _i in range(n_electrons - (n_qubits // 2))] + return sorted(indices) + + +def a_gate_indices(n_qubits, n_electrons): + """ + Obtain the qubit indices for a single layer of the primitive pattern of 'A' gates in the circuit ansatz """ - result = [] - result.append(gates.CNOT(qubit2, qubit1)) - result += [gates.RY(qubit2, 0.0), gates.RZ(qubit2, 0.0)] - result.append(gates.CNOT(qubit1, qubit2)) - result += [gates.RZ(qubit2, 0.0), gates.RY(qubit2, 0.0)] - result.append(gates.CNOT(qubit2, qubit1)) - return result + # 1. Apply X gates to qubits. Avoid placing gates on neighboring qubits + x_gates = x_gate_indices(n_qubits, n_electrons) + # 2. Apply 'first layer' of gates on all adjacent pairs of qubits on which either X*I or I*X has been applied. + first_layer = [(_i, _i + 1) for _i in x_gates if _i + 1 < n_qubits and _i + 1 not in x_gates] + first_layer += [(_i - 1, _i) for _i in x_gates if _i - 1 >= 0 and _i - 1 not in x_gates] + # 3a. Apply 'second layer' of gates on adjacent pairs of qubits. Each pair includes 1 qubit acted on in the previous + # step and a qubit free of gates. Continue placing gates on adjacent qubits until all neighboring qubits are connected + second_layer = [(_i, _i + 1) for _i in range(max(pair[1] for pair in first_layer), n_qubits - 1)] + second_layer += [(_i - 1, _i) for _i in range(min(pair[0] for pair in first_layer), 0, -1)] + # 3b. The first and second layers define a primitive pattern: + primitive_pattern = first_layer + second_layer + # 4. Repeat the primitive pattern until (n_qubits choose n_electrons) A gates are placed + n_gates_per_layer = len(primitive_pattern) + n_a_gates = n_qubits * (n_qubits - 1) // 2 + assert ( + n_a_gates % n_gates_per_layer == 0 + ), f"n_a_gates ({n_a_gates}) is not a multiple of n_gates_per_layer ({n_gates_per_layer})!" + return (n_a_gates // n_gates_per_layer) * primitive_pattern # Main function @@ -34,19 +85,10 @@ def symm_preserving_circuit(n_qubits, n_electrons): Qibo ``Circuit``: Circuit ansatz """ circuit = Circuit(n_qubits) - circuit.add(gates.X(2 * _i) for _i in range(n_electrons)) - - a_gate_qubits = [] # Generate the list of qubits pairs for adding A gates - - a_gate_qubits = [(0, 1)] - + circuit.add(gates.X(_i) for _i in x_gate_indices(n_qubits, n_electrons)) + # Generate the qubit pair indices for adding A gates + a_gate_qubits = a_gate_indices(n_qubits, n_electrons) a_gates = [a_gate(qubit1, qubit2) for qubit1, qubit2 in a_gate_qubits] - circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate) # Unpack the nested list - + # Each a_gate is a list of elementary gates, so a_gates is a nested list; need to unpack it + circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate) return circuit - - -n_qubits = 4 -n_electrons = 2 -circuit = symm_preserving_circuit(n_qubits, n_electrons) -print(circuit.draw()) From 5d2705a267362adee82f406c497565380e12233e Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Fri, 28 Jun 2024 06:53:00 +0000 Subject: [PATCH 04/32] Fix small bugs, should be working now Next: Write tests --- src/qibochem/ansatz/symmetry.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py index a753721..f85dc8e 100644 --- a/src/qibochem/ansatz/symmetry.py +++ b/src/qibochem/ansatz/symmetry.py @@ -2,6 +2,8 @@ Symmetry-Preserving circuit ansatz from Gard et al. Reference: https://doi.org/10.1038/s41534-019-0240-1 """ +from math import factorial + import numpy as np from qibo import Circuit, gates @@ -48,12 +50,10 @@ def x_gate_indices(n_qubits, n_electrons): return sorted(indices) -def a_gate_indices(n_qubits, n_electrons): +def a_gate_indices(n_qubits, n_electrons, x_gates): """ Obtain the qubit indices for a single layer of the primitive pattern of 'A' gates in the circuit ansatz """ - # 1. Apply X gates to qubits. Avoid placing gates on neighboring qubits - x_gates = x_gate_indices(n_qubits, n_electrons) # 2. Apply 'first layer' of gates on all adjacent pairs of qubits on which either X*I or I*X has been applied. first_layer = [(_i, _i + 1) for _i in x_gates if _i + 1 < n_qubits and _i + 1 not in x_gates] first_layer += [(_i - 1, _i) for _i in x_gates if _i - 1 >= 0 and _i - 1 not in x_gates] @@ -63,9 +63,11 @@ def a_gate_indices(n_qubits, n_electrons): second_layer += [(_i - 1, _i) for _i in range(min(pair[0] for pair in first_layer), 0, -1)] # 3b. The first and second layers define a primitive pattern: primitive_pattern = first_layer + second_layer + # Need to add any missing connections between neighbouring qubits + primitive_pattern += [pair for _i in range(n_qubits - 1) if (pair := (_i, _i + 1)) not in primitive_pattern] # 4. Repeat the primitive pattern until (n_qubits choose n_electrons) A gates are placed n_gates_per_layer = len(primitive_pattern) - n_a_gates = n_qubits * (n_qubits - 1) // 2 + n_a_gates = factorial(n_qubits) // (factorial(n_qubits - n_electrons) * factorial(n_electrons)) assert ( n_a_gates % n_gates_per_layer == 0 ), f"n_a_gates ({n_a_gates}) is not a multiple of n_gates_per_layer ({n_gates_per_layer})!" @@ -85,9 +87,10 @@ def symm_preserving_circuit(n_qubits, n_electrons): Qibo ``Circuit``: Circuit ansatz """ circuit = Circuit(n_qubits) - circuit.add(gates.X(_i) for _i in x_gate_indices(n_qubits, n_electrons)) + x_gates = x_gate_indices(n_qubits, n_electrons) + circuit.add(gates.X(_i) for _i in x_gates) # Generate the qubit pair indices for adding A gates - a_gate_qubits = a_gate_indices(n_qubits, n_electrons) + a_gate_qubits = a_gate_indices(n_qubits, n_electrons, x_gates) a_gates = [a_gate(qubit1, qubit2) for qubit1, qubit2 in a_gate_qubits] # Each a_gate is a list of elementary gates, so a_gates is a nested list; need to unpack it circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate) From ac3c2363bf7982e5fe788754c1376a6d80247f41 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Mon, 1 Jul 2024 07:14:38 +0000 Subject: [PATCH 05/32] Fix bugs and start adding tests --- src/qibochem/ansatz/__init__.py | 1 + src/qibochem/ansatz/symmetry.py | 13 +++++-- tests/test_symmetry_preserving_ansatz.py | 45 ++++++++++++++++++++++++ 3 files changed, 56 insertions(+), 3 deletions(-) create mode 100644 tests/test_symmetry_preserving_ansatz.py diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py index 02cee9c..8dbd59d 100644 --- a/src/qibochem/ansatz/__init__.py +++ b/src/qibochem/ansatz/__init__.py @@ -1,6 +1,7 @@ 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.symmetry import symm_preserving_circuit from qibochem.ansatz.ucc import ( generate_excitations, mp2_amplitude, diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py index f85dc8e..019ab65 100644 --- a/src/qibochem/ansatz/symmetry.py +++ b/src/qibochem/ansatz/symmetry.py @@ -38,8 +38,15 @@ def a_gate(qubit1, qubit2, theta=None, phi=None): phi = 0.0 # R(theta, phi) = R_z (phi + pi) R_y (theta + 0.5*pi) - r_and_cnot = [gates.RZ(qubit2, phi + np.pi), gates.RY(qubit2, theta + 0.5 * np.pi), gates.CNOT(qubit2, qubit1)] - return r_and_cnot[::-1] + [gates.CNOT(qubit1, qubit2)] + r_and_cnot + r_gate = [gates.RY(qubit2, theta + 0.5 * np.pi), gates.RZ(qubit2, phi + np.pi)] + r_gate_dagger = [_gate.dagger() for _gate in r_gate][::-1] + return ( + [gates.CNOT(qubit2, qubit1)] + + r_gate_dagger + + [gates.CNOT(qubit1, qubit2)] + + r_gate + + [gates.CNOT(qubit2, qubit1)] + ) def x_gate_indices(n_qubits, n_electrons): @@ -77,7 +84,7 @@ def a_gate_indices(n_qubits, n_electrons, x_gates): # Main function def symm_preserving_circuit(n_qubits, n_electrons): """ - Symmetry-preserving circuit ansatz + Symmetry-preserving circuit ansatz from Gard et al. Reference: https://doi.org/10.1038/s41534-019-0240-1 Args: n_qubits: Number of qubits in the quantum circuit diff --git a/tests/test_symmetry_preserving_ansatz.py b/tests/test_symmetry_preserving_ansatz.py new file mode 100644 index 0000000..c5026aa --- /dev/null +++ b/tests/test_symmetry_preserving_ansatz.py @@ -0,0 +1,45 @@ +import numpy as np +import pytest +from qibo import Circuit, gates + +from qibochem.ansatz import symm_preserving_circuit +from qibochem.ansatz.symmetry import a_gate, a_gate_indices, x_gate_indices + + +@pytest.mark.parametrize( + "theta,phi,expected", + [ + (0.0, 0.0, np.array([0.0, 0.0, -1.0, 00])), + (0.5 * np.pi, 0.0, np.array([0.0, 1.0, 0.0, 00])), + (0.5 * np.pi, np.pi, np.array([0.0, -1.0, 0.0, 00])), + ], +) +def test_a_gate(theta, phi, expected): + circuit = Circuit(2) + circuit.add(gates.X(0)) + a_gates = a_gate(0, 1, theta=theta, phi=phi) + circuit.add(a_gates) + + result = circuit(nshots=1) + state_ket = result.state() + assert np.allclose(state_ket, expected) + + +def test_x_gate_indices(): + # x_gate_indices(n_qubits, n_electrons) + pass + + +def test_a_gate_indices(): + # a_gate_indices(n_qubits, n_electrons, x_gates + pass + + +def test_symm_preserving_circuit(): + # symm_preserving_circuit(n_qubits, n_electrons + pass + + +def test_vqe_symm_preserving_circuit(): + # Maybe don't do...? + pass From 1314849142e972d1f4566204082d21848c6f79ee Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Mon, 1 Jul 2024 08:21:41 +0000 Subject: [PATCH 06/32] Add more tests for symmetry ansatz --- src/qibochem/ansatz/symmetry.py | 1 + tests/test_symmetry_preserving_ansatz.py | 66 ++++++++++++++++++------ 2 files changed, 52 insertions(+), 15 deletions(-) diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py index 019ab65..170df7b 100644 --- a/src/qibochem/ansatz/symmetry.py +++ b/src/qibochem/ansatz/symmetry.py @@ -61,6 +61,7 @@ def a_gate_indices(n_qubits, n_electrons, x_gates): """ Obtain the qubit indices for a single layer of the primitive pattern of 'A' gates in the circuit ansatz """ + assert len(x_gates) == n_electrons, f"n_electrons ({n_electrons}) != Number of X gates given! ({x_gates})" # 2. Apply 'first layer' of gates on all adjacent pairs of qubits on which either X*I or I*X has been applied. first_layer = [(_i, _i + 1) for _i in x_gates if _i + 1 < n_qubits and _i + 1 not in x_gates] first_layer += [(_i - 1, _i) for _i in x_gates if _i - 1 >= 0 and _i - 1 not in x_gates] diff --git a/tests/test_symmetry_preserving_ansatz.py b/tests/test_symmetry_preserving_ansatz.py index c5026aa..72a7fa7 100644 --- a/tests/test_symmetry_preserving_ansatz.py +++ b/tests/test_symmetry_preserving_ansatz.py @@ -1,15 +1,23 @@ +""" +Tests for the symmetry preserving ansatz from Gard et al. (DOI: https://doi.org/10.1038/s41534-019-0240-1) +""" + import numpy as np import pytest from qibo import Circuit, gates -from qibochem.ansatz import symm_preserving_circuit -from qibochem.ansatz.symmetry import a_gate, a_gate_indices, x_gate_indices +from qibochem.ansatz.symmetry import ( + a_gate, + a_gate_indices, + symm_preserving_circuit, + x_gate_indices, +) @pytest.mark.parametrize( "theta,phi,expected", [ - (0.0, 0.0, np.array([0.0, 0.0, -1.0, 00])), + (None, None, np.array([0.0, 0.0, -1.0, 00])), (0.5 * np.pi, 0.0, np.array([0.0, 1.0, 0.0, 00])), (0.5 * np.pi, np.pi, np.array([0.0, -1.0, 0.0, 00])), ], @@ -25,21 +33,49 @@ def test_a_gate(theta, phi, expected): assert np.allclose(state_ket, expected) -def test_x_gate_indices(): - # x_gate_indices(n_qubits, n_electrons) - pass +@pytest.mark.parametrize( + "n_qubits,n_electrons,expected", + [ + (4, 2, [0, 2]), + (4, 3, [0, 1, 2]), + (4, 4, [0, 1, 2, 3]), + ], +) +def test_x_gate_indices(n_qubits, n_electrons, expected): + test = x_gate_indices(n_qubits, n_electrons) + assert test == expected -def test_a_gate_indices(): - # a_gate_indices(n_qubits, n_electrons, x_gates - pass +@pytest.mark.parametrize( + "n_qubits,n_electrons,x_gates,expected", + [ + (4, 2, [0, 2], 2 * [(0, 1), (2, 3), (1, 2)]), + (6, 4, [0, 1, 2, 4], 3 * [(2, 3), (4, 5), (3, 4), (1, 2), (0, 1)]), + ], +) +def test_a_gate_indices(n_qubits, n_electrons, x_gates, expected): + test = a_gate_indices(n_qubits, n_electrons, x_gates) + assert test == expected -def test_symm_preserving_circuit(): - # symm_preserving_circuit(n_qubits, n_electrons - pass +@pytest.mark.parametrize( + "n_qubits,n_electrons", + [ + (4, 2), + (6, 4), + ], +) +def test_symm_preserving_circuit(n_qubits, n_electrons): + control_circuit = Circuit(n_qubits) + x_gates = x_gate_indices(n_qubits, n_electrons) + control_circuit.add(gates.X(_i) for _i in x_gates) + a_gate_qubits = a_gate_indices(n_qubits, n_electrons, x_gates) + a_gates = [a_gate(qubit1, qubit2) for qubit1, qubit2 in a_gate_qubits] + control_circuit.add(_gates for _a_gate in a_gates for _gates in _a_gate) + test_circuit = symm_preserving_circuit(n_qubits, n_electrons) -def test_vqe_symm_preserving_circuit(): - # Maybe don't do...? - pass + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(list(control_circuit.queue), list(test_circuit.queue)) + ) From 9b606104be382de9ba245e287f9806434a6df909 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Mon, 29 Jul 2024 08:36:48 +0000 Subject: [PATCH 07/32] Remove old draft file Ansatz for paper 2 should be done, no need to keep the initial draft anymore --- examples/paper2.py | 58 ---------------------------------------------- 1 file changed, 58 deletions(-) delete mode 100644 examples/paper2.py diff --git a/examples/paper2.py b/examples/paper2.py deleted file mode 100644 index 3159d6b..0000000 --- a/examples/paper2.py +++ /dev/null @@ -1,58 +0,0 @@ -""" -A rough Python script attempting to implement the VQE circuit ansatz proposed in Gard et al. using Qibo. - -Reference paper: https://doi.org/10.1038/s41534-019-0240-1 - -Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Caleb Seow from -Eunoia Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for -Junior College students. -""" - -import numpy as np -from qibo import Circuit, gates - -from qibochem.driver.molecule import Molecule - -# Define molecule -mol = Molecule(xyz_file="h2.xyz") -mol.run_pyscf() - - -n = 4 # no of qubits -m = 2 # no of electrons -theta = 0.1 -phi = 0.1 - - -def special_R(a, theta, phi): - A = Circuit(n) - A.add(gates.RY(a, theta + np.pi / 2)) - A.add(gates.RZ(a, phi + np.pi)) - return A - - -print(special_R(1, 0.1, 0.1).invert().draw()) -print(special_R(1, 0.1, 0.1).draw()) - - -def A_gate(a): - A = Circuit(n) - A.add(gates.CNOT(a + 1, a)) - A += special_R(a + 1, theta, phi).invert() - A.add(gates.CNOT(a, a + 1)) - A += special_R(a + 1, theta, phi) - A.add(gates.CNOT(a + 1, a)) - return A - - -print(A_gate(0).draw()) - -Gate2 = Circuit(n) -Gate2.add(gates.X(1)) -Gate2.add(gates.X(2)) -for i in range(2): - Gate2 += A_gate(0) - Gate2 += A_gate(2) - Gate2 += A_gate(1) - -print(Gate2.draw()) From 65fe973b0c61fba73f1f1aaf5bb2551952751756 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Mon, 29 Jul 2024 09:38:52 +0000 Subject: [PATCH 08/32] Draft code for remaining circuit ansatz --- src/qibochem/ansatz/universal.py | 83 ++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 src/qibochem/ansatz/universal.py diff --git a/src/qibochem/ansatz/universal.py b/src/qibochem/ansatz/universal.py new file mode 100644 index 0000000..0573b4b --- /dev/null +++ b/src/qibochem/ansatz/universal.py @@ -0,0 +1,83 @@ +""" +Module documentation +""" + +from qibo import Circuit, gates + +# Helper functions + + +def double_excitation_gate(excitation, theta=0.0): + """ + Decomposition of a double excitation gate into single qubit rotations and CNOTs + + Args: + n_qubits (int): Number of qubits in the circuit + excitation: Iterable of orbitals involved in the excitation; must have an even number of elements + theta (float): Rotation angle. Default: 0.0 + + Returns: + (list): List of gates representing the decomposition of the Givens' double excitation gate + """ + sorted_orbitals = sorted(excitation) + # Check size of orbitals input + assert len(sorted_orbitals) % 2 == 0, f"{excitation} must have an even number of items" + + if theta is None: + theta = 0.0 + + result = [] + result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3])) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[2])) + result.append(gates.H(sorted_orbitals[0])) + result.append(gates.H(sorted_orbitals[3])) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1])) + result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3])) + result.append(gates.RY(sorted_orbitals[0], -0.125 * theta)) + result.append(gates.RY(sorted_orbitals[1], 0.125 * theta)) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[3])) + result.append(gates.H(sorted_orbitals[3])) + result.append(gates.CNOT(sorted_orbitals[3], sorted_orbitals[1])) + result.append(gates.RY(sorted_orbitals[0], -0.125 * theta)) + result.append(gates.RY(sorted_orbitals[1], 0.125 * theta)) + result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[1])) + result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[0])) + result.append(gates.RY(sorted_orbitals[0], 0.125 * theta)) + result.append(gates.RY(sorted_orbitals[1], -0.125 * theta)) + result.append(gates.CNOT(sorted_orbitals[3], sorted_orbitals[1])) + result.append(gates.H(sorted_orbitals[3])) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[3])) + result.append(gates.RY(sorted_orbitals[0], 0.125 * theta)) + result.append(gates.RY(sorted_orbitals[1], -0.125 * theta)) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1])) + result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[0])) + result.append(gates.H(sorted_orbitals[0])) + result.append(gates.H(sorted_orbitals[3])) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[2])) + result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3])) + return result + + +# Main function +def universal_circuit(n_qubits, excitation, theta=0.0): + """ + Universal circuit ansatz from Arrazola et al. Reference: https://doi.org/10.22331/q-2022-06-20-742 + + Args: + n_qubits: Number of qubits in the quantum circuit + n_electrons: Number of electrons in the molecular system + + Returns: + Qibo ``Circuit``: Circuit ansatz + """ + circuit = Circuit(n_qubits) + if len(excitation) == 2: + circuit.add(gates.GIVENS(excitation[0], excitation[1], theta)) + elif len(excitation) == 4: + circuit.add(double_excitation_gate(excitation, theta)) + return circuit + + +def universal_ansatz(n_qubits, excitation, n_electrons): + """TODO: Same implementation as UCC?""" + pass From d5449e71092dc0250d8085ba93c1de749a515d5f Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 31 Jul 2024 06:52:45 +0000 Subject: [PATCH 09/32] Finish code for circuit with only one excitation i.e. return a circuit with either a single or double excitation represented as a Givens rotation --- src/qibochem/ansatz/universal.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/src/qibochem/ansatz/universal.py b/src/qibochem/ansatz/universal.py index 0573b4b..9ad930e 100644 --- a/src/qibochem/ansatz/universal.py +++ b/src/qibochem/ansatz/universal.py @@ -7,22 +7,17 @@ # Helper functions -def double_excitation_gate(excitation, theta=0.0): +def double_excitation_gate(sorted_orbitals, theta=0.0): """ - Decomposition of a double excitation gate into single qubit rotations and CNOTs + Decomposition of a Givens double excitation gate into single qubit rotations and CNOTs Args: - n_qubits (int): Number of qubits in the circuit - excitation: Iterable of orbitals involved in the excitation; must have an even number of elements + sorted_orbitals (list): Sorted list of orbitals involved in the excitation theta (float): Rotation angle. Default: 0.0 Returns: (list): List of gates representing the decomposition of the Givens' double excitation gate """ - sorted_orbitals = sorted(excitation) - # Check size of orbitals input - assert len(sorted_orbitals) % 2 == 0, f"{excitation} must have an even number of items" - if theta is None: theta = 0.0 @@ -59,22 +54,28 @@ def double_excitation_gate(excitation, theta=0.0): # Main function -def universal_circuit(n_qubits, excitation, theta=0.0): +def givens_excitation_circuit(n_qubits, excitation, theta=0.0): """ - Universal circuit ansatz from Arrazola et al. Reference: https://doi.org/10.22331/q-2022-06-20-742 + Circuit ansatz for one Givens rotation excitation from Arrazola et al. Reference: + https://doi.org/10.22331/q-2022-06-20-742 Args: n_qubits: Number of qubits in the quantum circuit - n_electrons: Number of electrons in the molecular system + excitation: Iterable of orbitals involved in the excitation; must have an even number of elements + E.g. ``[0, 1, 2, 3]`` represents the excitation of electrons in orbitals ``(0, 1)`` to ``(2, 3)`` Returns: Qibo ``Circuit``: Circuit ansatz """ + sorted_orbitals = sorted(excitation) + # Check size of orbitals input + assert len(sorted_orbitals) % 2 == 0, f"{excitation} must have an even number of items" + circuit = Circuit(n_qubits) if len(excitation) == 2: circuit.add(gates.GIVENS(excitation[0], excitation[1], theta)) elif len(excitation) == 4: - circuit.add(double_excitation_gate(excitation, theta)) + circuit.add(double_excitation_gate(sorted_orbitals, theta)) return circuit From 4c88ffd8ed66cdc2f283786312275ffe8582c093 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 31 Jul 2024 06:57:10 +0000 Subject: [PATCH 10/32] Add error handling if excitation input not valid --- src/qibochem/ansatz/universal.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/qibochem/ansatz/universal.py b/src/qibochem/ansatz/universal.py index 9ad930e..a8e3bef 100644 --- a/src/qibochem/ansatz/universal.py +++ b/src/qibochem/ansatz/universal.py @@ -76,6 +76,8 @@ def givens_excitation_circuit(n_qubits, excitation, theta=0.0): circuit.add(gates.GIVENS(excitation[0], excitation[1], theta)) elif len(excitation) == 4: circuit.add(double_excitation_gate(sorted_orbitals, theta)) + else: + raise NotImplementedError("Can only handle single and double excitations!") return circuit From 0f69ccf1007d38b50d925da182f29cc748511e72 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 31 Jul 2024 06:58:01 +0000 Subject: [PATCH 11/32] Rename file; 'universal.py' is too vague --- src/qibochem/ansatz/{universal.py => givens_excitation.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/qibochem/ansatz/{universal.py => givens_excitation.py} (100%) diff --git a/src/qibochem/ansatz/universal.py b/src/qibochem/ansatz/givens_excitation.py similarity index 100% rename from src/qibochem/ansatz/universal.py rename to src/qibochem/ansatz/givens_excitation.py From c181fb12da3c51d7e86454661fc4837a10c133f3 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 31 Jul 2024 07:41:12 +0000 Subject: [PATCH 12/32] Move out some utility functions to util.py --- src/qibochem/ansatz/__init__.py | 11 +-- src/qibochem/ansatz/ucc.py | 1 + src/qibochem/ansatz/util.py | 154 ++++++++++++++++++++++++++++++++ 3 files changed, 157 insertions(+), 9 deletions(-) create mode 100644 src/qibochem/ansatz/util.py diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py index 8dbd59d..e0b77f5 100644 --- a/src/qibochem/ansatz/__init__.py +++ b/src/qibochem/ansatz/__init__.py @@ -2,12 +2,5 @@ from qibochem.ansatz.hardware_efficient import he_circuit from qibochem.ansatz.hf_reference import hf_circuit from qibochem.ansatz.symmetry import symm_preserving_circuit -from qibochem.ansatz.ucc import ( - 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' +from qibochem.ansatz.ucc import ucc_ansatz, ucc_circuit +from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations diff --git a/src/qibochem/ansatz/ucc.py b/src/qibochem/ansatz/ucc.py index 4fde184..2597bc4 100644 --- a/src/qibochem/ansatz/ucc.py +++ b/src/qibochem/ansatz/ucc.py @@ -7,6 +7,7 @@ from qibo import Circuit, gates from qibochem.ansatz.hf_reference import hf_circuit +from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations def expi_pauli(n_qubits, pauli_string, theta): diff --git a/src/qibochem/ansatz/util.py b/src/qibochem/ansatz/util.py new file mode 100644 index 0000000..a3521af --- /dev/null +++ b/src/qibochem/ansatz/util.py @@ -0,0 +1,154 @@ +""" +Utility functions that can be used by different ansatzes +""" + + +def mp2_amplitude(excitation, orbital_energies, tei): + r""" + Calculate the MP2 guess amplitude for a single UCC circuit: 0.0 for a single excitation. + for a double excitation (In SO basis): :math:`t_{ij}^{ab} = (g_{ijab} - g_{ijba}) / (e_i + e_j - e_a - e_b)` + + Args: + excitation: Iterable of spin-orbitals representing a excitation. Must have either 2 or 4 elements exactly, + representing a single or double excitation respectively. + orbital_energies: eigenvalues of the Fock operator, i.e. orbital energies + tei: Two-electron integrals in MO basis and second quantization notation + + Returns: + MP2 guess amplitude (float) + """ + # Checks validity of excitation argument + assert len(excitation) % 2 == 0 and len(excitation) // 2 <= 2, f"{excitation} must have either 2 or 4 elements" + # If single excitation, can just return 0.0 directly + if len(excitation) == 2: + return 0.0 + + # Convert orbital indices to be in MO basis + mo_orbitals = [orbital // 2 for orbital in excitation] + # Numerator: g_ijab - g_ijba + g_ijab = ( + tei[tuple(mo_orbitals)] # Can index directly using the MO TEIs + if (excitation[0] + excitation[3]) % 2 == 0 and (excitation[1] + excitation[2]) % 2 == 0 + else 0.0 + ) + g_ijba = ( + tei[tuple(mo_orbitals[:2] + mo_orbitals[2:][::-1])] # Reverse last two terms + if (excitation[0] + excitation[2]) % 2 == 0 and (excitation[1] + excitation[3]) % 2 == 0 + else 0.0 + ) + numerator = g_ijab - g_ijba + # Denominator is directly from the orbital energies + denominator = sum(orbital_energies[mo_orbitals[:2]]) - sum(orbital_energies[mo_orbitals[2:]]) + return numerator / denominator + + +def generate_excitations(order, excite_from, excite_to, conserve_spin=True): + """ + Generate all possible excitations between a list of occupied and virtual orbitals + + Args: + order: Order of excitations, i.e. 1 == single, 2 == double + excite_from: Iterable of integers + excite_to: Iterable of integers + conserve_spin: ensure that the total electronic spin is conserved + + Return: + List of lists, e.g. [[0, 1]] + """ + # If order of excitation > either number of electrons/orbitals, return list of empty list + if order > min(len(excite_from), len(excite_to)): + return [[]] + + # Generate all possible excitations first + 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) + ] + # Filter out the excitations if conserve_spin set + if conserve_spin: + # Not sure if this filtering is exhaustive; might not remove some redundant excitations? + all_excitations = [ + _ex + for _ex in all_excitations + if sum(_ex) % 2 == 0 and (sum(_i % 2 for _i in _ex[:order]) == sum(_i % 2 for _i in _ex[order:])) + ] + return all_excitations + + +def sort_excitations(excitations): + """ + Sorts a list of excitations according to some common-sense and empirical rules (see below). + The order of the excitations must be the same throughout. + + Sorting order: + 1. (For double excitations only) All paired excitations between the same MOs first + 2. Pair up excitations between the same MOs, e.g. (0, 2) and (1, 3) + 3. Then count upwards from smallest MO + + Args: + excitations: List of iterables, each representing an excitation; e.g. [[1, 5], [0, 4]] + + Returns: + List of excitations after sorting + """ + # Check that all excitations are of the same order and <= 2 + order = len(excitations[0]) // 2 + if order > 2: + raise NotImplementedError("Can only handle single and double excitations!") + assert all(len(_ex) // 2 == order for _ex in excitations), "Cannot handle excitations of different orders!" + + # Define variables for the while loop + copy_excitations = [list(_ex) for _ex in excitations] + result = [] + prev = [] + + # No idea how I came up with this, but it seems to work for double excitations + 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: + if not prev: + # Take out all pair excitations first + pair_excitations = [ + _ex + for _ex in copy_excitations + # Indices of the electrons/holes must be consecutive numbers + if sum(abs(_ex[2 * _i + 1] // 2 - _ex[2 * _i] // 2) for _i in range(0, order)) == 0 + ] + while pair_excitations: + pair_excitations = sorted(pair_excitations) + ex_to_remove = pair_excitations.pop(0) + if ex_to_remove in copy_excitations: + # 'Move' the first pair excitation from copy_excitations to result + index = copy_excitations.index(ex_to_remove) + result.append(copy_excitations.pop(index)) + + # No more pair excitations, only remaining excitations should have >=3 MOs involved + # Sort the remaining excitations + copy_excitations = sorted(copy_excitations, key=sorting_fn if order != 1 else None) + else: + # Check to see for excitations involving the same MOs as prev + _from = prev[:order] + _to = prev[order:] + # Get all possible excitations involving the same MOs + new_from = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _from] + new_to = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _to] + same_mo_ex = [sorted(list(_f) + list(_t)) for _f in (_from, new_from) for _t in (_to, new_to)] + # Remove the excitations with the same MOs from copy_excitations + while same_mo_ex: + same_mo_ex = sorted(same_mo_ex) + ex_to_remove = same_mo_ex.pop(0) + if ex_to_remove in copy_excitations: + # 'Move' the first entry of same_mo_index from copy_excitations to result + index = copy_excitations.index(ex_to_remove) + result.append(copy_excitations.pop(index)) + prev = None + continue + if copy_excitations: + # Remove the first entry from the sorted list of remaining excitations and add it to result + prev = copy_excitations.pop(0) + result.append(prev) + return result From 2098761a9d2318767df55141227ca1929b27caf2 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 31 Jul 2024 07:53:11 +0000 Subject: [PATCH 13/32] Move utility function tests over to new .py file --- tests/test_ansatz_util.py | 53 ++++++++++++++++++++++++++++++++++++ tests/test_ucc.py | 56 +++------------------------------------ 2 files changed, 56 insertions(+), 53 deletions(-) create mode 100644 tests/test_ansatz_util.py diff --git a/tests/test_ansatz_util.py b/tests/test_ansatz_util.py new file mode 100644 index 0000000..9a81548 --- /dev/null +++ b/tests/test_ansatz_util.py @@ -0,0 +1,53 @@ +""" +Tests for ansatz utility functions (util.py) +""" + +import numpy as np +import pytest + +from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations +from qibochem.driver import Molecule + + +@pytest.mark.parametrize( + "order,excite_from,excite_to,expected", + [ + (1, [0, 1], [2, 3], [[0, 2], [1, 3]]), + (2, [2, 3], [4, 5], [[2, 3, 4, 5]]), + (3, [0, 1], [2, 3], [[]]), + ], +) +def test_generate_excitations(order, excite_from, excite_to, expected): + """Test generation of all possible excitations between two lists of orbitals""" + test = generate_excitations(order, excite_from, excite_to) + assert test == expected + + +@pytest.mark.parametrize( + "order,excite_from,excite_to,expected", + [ + (1, [0, 1], [2, 3, 4, 5], [[0, 2], [1, 3], [0, 4], [1, 5]]), + (2, [0, 1], [2, 3, 4, 5], [[0, 1, 2, 3], [0, 1, 4, 5], [0, 1, 2, 5], [0, 1, 3, 4]]), + ], +) +def test_sort_excitations(order, excite_from, excite_to, expected): + test = sort_excitations(generate_excitations(order, excite_from, excite_to)) + assert test == expected + + +def test_sort_excitations_triples(): + with pytest.raises(NotImplementedError): + sort_excitations([[1, 2, 3, 4, 5, 6]]) + + +def test_mp2_amplitude_singles(): + assert mp2_amplitude([0, 2], np.random.rand(4), np.random.rand(4, 4)) == 0.0 + + +def test_mp2_amplitude_doubles(): + h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) + h2.run_pyscf() + l = mp2_amplitude([0, 1, 2, 3], h2.eps, h2.tei) + ref_l = 0.06834019757197053 + + assert np.isclose(l, ref_l) diff --git a/tests/test_ucc.py b/tests/test_ucc.py index d6d6cc7..5b7ef11 100644 --- a/tests/test_ucc.py +++ b/tests/test_ucc.py @@ -1,5 +1,5 @@ """ -Tests for the UCC ansatz and related functions +Tests for the UCC ansatz """ from functools import reduce @@ -10,61 +10,11 @@ from qibo.hamiltonians import SymbolicHamiltonian from qibochem.ansatz import hf_circuit -from qibochem.ansatz.ucc import ( - expi_pauli, - generate_excitations, - mp2_amplitude, - sort_excitations, - ucc_ansatz, - ucc_circuit, -) +from qibochem.ansatz.ucc import expi_pauli, ucc_ansatz, ucc_circuit +from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations from qibochem.driver import Molecule -@pytest.mark.parametrize( - "order,excite_from,excite_to,expected", - [ - (1, [0, 1], [2, 3], [[0, 2], [1, 3]]), - (2, [2, 3], [4, 5], [[2, 3, 4, 5]]), - (3, [0, 1], [2, 3], [[]]), - ], -) -def test_generate_excitations(order, excite_from, excite_to, expected): - """Test generation of all possible excitations between two lists of orbitals""" - test = generate_excitations(order, excite_from, excite_to) - assert test == expected - - -@pytest.mark.parametrize( - "order,excite_from,excite_to,expected", - [ - (1, [0, 1], [2, 3, 4, 5], [[0, 2], [1, 3], [0, 4], [1, 5]]), - (2, [0, 1], [2, 3, 4, 5], [[0, 1, 2, 3], [0, 1, 4, 5], [0, 1, 2, 5], [0, 1, 3, 4]]), - ], -) -def test_sort_excitations(order, excite_from, excite_to, expected): - test = sort_excitations(generate_excitations(order, excite_from, excite_to)) - assert test == expected - - -def test_sort_excitations_triples(): - with pytest.raises(NotImplementedError): - sort_excitations([[1, 2, 3, 4, 5, 6]]) - - -def test_mp2_amplitude_singles(): - assert mp2_amplitude([0, 2], np.random.rand(4), np.random.rand(4, 4)) == 0.0 - - -def test_mp2_amplitude_doubles(): - h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) - h2.run_pyscf() - l = mp2_amplitude([0, 1, 2, 3], h2.eps, h2.tei) - ref_l = 0.06834019757197053 - - assert np.isclose(l, ref_l) - - @pytest.mark.parametrize( "pauli_string", [ From f5d30206c8a3941e2bbc2c979ced5487652107bf Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 31 Jul 2024 08:02:52 +0000 Subject: [PATCH 14/32] Remove utility functions from ucc.py --- src/qibochem/ansatz/ucc.py | 154 ------------------------------------- 1 file changed, 154 deletions(-) diff --git a/src/qibochem/ansatz/ucc.py b/src/qibochem/ansatz/ucc.py index 2597bc4..99a9b78 100644 --- a/src/qibochem/ansatz/ucc.py +++ b/src/qibochem/ansatz/ucc.py @@ -109,160 +109,6 @@ def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map return circuit -# Utility functions - - -def mp2_amplitude(excitation, orbital_energies, tei): - r""" - Calculate the MP2 guess amplitude for a single UCC circuit: 0.0 for a single excitation. - for a double excitation (In SO basis): :math:`t_{ij}^{ab} = (g_{ijab} - g_{ijba}) / (e_i + e_j - e_a - e_b)` - - Args: - excitation: Iterable of spin-orbitals representing a excitation. Must have either 2 or 4 elements exactly, - representing a single or double excitation respectively. - orbital_energies: eigenvalues of the Fock operator, i.e. orbital energies - tei: Two-electron integrals in MO basis and second quantization notation - - Returns: - MP2 guess amplitude (float) - """ - # Checks validity of excitation argument - assert len(excitation) % 2 == 0 and len(excitation) // 2 <= 2, f"{excitation} must have either 2 or 4 elements" - # If single excitation, can just return 0.0 directly - if len(excitation) == 2: - return 0.0 - - # Convert orbital indices to be in MO basis - mo_orbitals = [orbital // 2 for orbital in excitation] - # Numerator: g_ijab - g_ijba - g_ijab = ( - tei[tuple(mo_orbitals)] # Can index directly using the MO TEIs - if (excitation[0] + excitation[3]) % 2 == 0 and (excitation[1] + excitation[2]) % 2 == 0 - else 0.0 - ) - g_ijba = ( - tei[tuple(mo_orbitals[:2] + mo_orbitals[2:][::-1])] # Reverse last two terms - if (excitation[0] + excitation[2]) % 2 == 0 and (excitation[1] + excitation[3]) % 2 == 0 - else 0.0 - ) - numerator = g_ijab - g_ijba - # Denominator is directly from the orbital energies - denominator = sum(orbital_energies[mo_orbitals[:2]]) - sum(orbital_energies[mo_orbitals[2:]]) - return numerator / denominator - - -def generate_excitations(order, excite_from, excite_to, conserve_spin=True): - """ - Generate all possible excitations between a list of occupied and virtual orbitals - - Args: - order: Order of excitations, i.e. 1 == single, 2 == double - excite_from: Iterable of integers - excite_to: Iterable of integers - conserve_spin: ensure that the total electronic spin is conserved - - Return: - List of lists, e.g. [[0, 1]] - """ - # If order of excitation > either number of electrons/orbitals, return list of empty list - if order > min(len(excite_from), len(excite_to)): - return [[]] - - # Generate all possible excitations first - 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) - ] - # Filter out the excitations if conserve_spin set - if conserve_spin: - # Not sure if this filtering is exhaustive; might not remove some redundant excitations? - all_excitations = [ - _ex - for _ex in all_excitations - if sum(_ex) % 2 == 0 and (sum(_i % 2 for _i in _ex[:order]) == sum(_i % 2 for _i in _ex[order:])) - ] - return all_excitations - - -def sort_excitations(excitations): - """ - Sorts a list of excitations according to some common-sense and empirical rules (see below). - The order of the excitations must be the same throughout. - - Sorting order: - 1. (For double excitations only) All paired excitations between the same MOs first - 2. Pair up excitations between the same MOs, e.g. (0, 2) and (1, 3) - 3. Then count upwards from smallest MO - - Args: - excitations: List of iterables, each representing an excitation; e.g. [[1, 5], [0, 4]] - - Returns: - List of excitations after sorting - """ - # Check that all excitations are of the same order and <= 2 - order = len(excitations[0]) // 2 - if order > 2: - raise NotImplementedError("Can only handle single and double excitations!") - assert all(len(_ex) // 2 == order for _ex in excitations), "Cannot handle excitations of different orders!" - - # Define variables for the while loop - copy_excitations = [list(_ex) for _ex in excitations] - result = [] - prev = [] - - # No idea how I came up with this, but it seems to work for double excitations - 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: - if not prev: - # Take out all pair excitations first - pair_excitations = [ - _ex - for _ex in copy_excitations - # Indices of the electrons/holes must be consecutive numbers - if sum(abs(_ex[2 * _i + 1] // 2 - _ex[2 * _i] // 2) for _i in range(0, order)) == 0 - ] - while pair_excitations: - pair_excitations = sorted(pair_excitations) - ex_to_remove = pair_excitations.pop(0) - if ex_to_remove in copy_excitations: - # 'Move' the first pair excitation from copy_excitations to result - index = copy_excitations.index(ex_to_remove) - result.append(copy_excitations.pop(index)) - - # No more pair excitations, only remaining excitations should have >=3 MOs involved - # Sort the remaining excitations - copy_excitations = sorted(copy_excitations, key=sorting_fn if order != 1 else None) - else: - # Check to see for excitations involving the same MOs as prev - _from = prev[:order] - _to = prev[order:] - # Get all possible excitations involving the same MOs - new_from = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _from] - new_to = [_i + 1 if _i % 2 == 0 else _i - 1 for _i in _to] - same_mo_ex = [sorted(list(_f) + list(_t)) for _f in (_from, new_from) for _t in (_to, new_to)] - # Remove the excitations with the same MOs from copy_excitations - while same_mo_ex: - same_mo_ex = sorted(same_mo_ex) - ex_to_remove = same_mo_ex.pop(0) - if ex_to_remove in copy_excitations: - # 'Move' the first entry of same_mo_index from copy_excitations to result - index = copy_excitations.index(ex_to_remove) - result.append(copy_excitations.pop(index)) - prev = None - continue - if copy_excitations: - # Remove the first entry from the sorted list of remaining excitations and add it to result - prev = copy_excitations.pop(0) - result.append(prev) - return result - - def ucc_ansatz( molecule, excitation_level=None, From 5b1b7258e41a006226890fd688cf52ca9f610075 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 03:08:29 +0000 Subject: [PATCH 15/32] Add tests for Givens excitation circuit --- src/qibochem/ansatz/givens_excitation.py | 2 +- tests/test_givens_excitation.py | 73 ++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 1 deletion(-) create mode 100644 tests/test_givens_excitation.py diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py index a8e3bef..c0192dc 100644 --- a/src/qibochem/ansatz/givens_excitation.py +++ b/src/qibochem/ansatz/givens_excitation.py @@ -81,6 +81,6 @@ def givens_excitation_circuit(n_qubits, excitation, theta=0.0): return circuit -def universal_ansatz(n_qubits, excitation, n_electrons): +def givens_excitation_ansatz(n_qubits, excitation, n_electrons): """TODO: Same implementation as UCC?""" pass diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py new file mode 100644 index 0000000..c1ee053 --- /dev/null +++ b/tests/test_givens_excitation.py @@ -0,0 +1,73 @@ +""" +Tests for the symmetry preserving ansatz from Gard et al. (DOI: https://doi.org/10.1038/s41534-019-0240-1) +""" + +import numpy as np +import pytest +from qibo import Circuit, gates + +from qibochem.ansatz.givens_excitation import ( + double_excitation_gate, + givens_excitation_ansatz, + givens_excitation_circuit, +) + + +def test_double_excitation_gate(): + # Hardcoded test + theta = 0.0 + + control_gates = [ + gates.CNOT(2, 3), + gates.CNOT(0, 2), + gates.H(0), + gates.H(3), + gates.CNOT(0, 1), + gates.CNOT(2, 3), + gates.RY(0, -theta), + gates.RY(1, theta), + gates.CNOT(0, 3), + gates.H(3), + gates.CNOT(3, 1), + gates.RY(0, -theta), + gates.RY(1, theta), + gates.CNOT(2, 1), + gates.CNOT(2, 0), + gates.RY(0, theta), + gates.RY(1, -theta), + gates.CNOT(3, 1), + gates.H(3), + gates.CNOT(0, 3), + gates.RY(0, theta), + gates.RY(1, -theta), + gates.CNOT(0, 1), + gates.CNOT(2, 0), + gates.H(0), + gates.H(3), + gates.CNOT(0, 2), + gates.CNOT(2, 3), + ] + + # Test circuit + test_circuit = givens_excitation_circuit(4, [0, 1, 2, 3]) + + # Check gates are correct + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(control_gates, list(test_circuit.queue)) + ) + + +@pytest.mark.parametrize( + "excitation,expected", + [ + ([0, 2], [gates.GIVENS(0, 2, 0.0)]), + ([0, 1, 2, 3], double_excitation_gate([0, 1, 2, 3])), + ], +) +def test_givens_excitation_circuit(excitation, expected): + test_circuit = givens_excitation_circuit(4, excitation) + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(expected, list(test_circuit.queue)) + ) From 1038d308d2a181993607c756359edc18a5c8d6f2 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 03:18:03 +0000 Subject: [PATCH 16/32] Fix code coverage? --- tests/test_givens_excitation.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py index c1ee053..682592d 100644 --- a/tests/test_givens_excitation.py +++ b/tests/test_givens_excitation.py @@ -47,14 +47,12 @@ def test_double_excitation_gate(): gates.CNOT(0, 2), gates.CNOT(2, 3), ] - - # Test circuit - test_circuit = givens_excitation_circuit(4, [0, 1, 2, 3]) + test_list = double_excitation_gate([0, 1, 2, 3]) # Check gates are correct assert all( control.name == test.name and control.target_qubits == test.target_qubits - for control, test in zip(control_gates, list(test_circuit.queue)) + for control, test in zip(control_gates, test_list) ) From c3847cf7d8ca9bebf67b67f2bf1217a96d1a3022 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 03:40:08 +0000 Subject: [PATCH 17/32] Fill in remaining code coverage --- tests/test_givens_excitation.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py index 682592d..769fa0d 100644 --- a/tests/test_givens_excitation.py +++ b/tests/test_givens_excitation.py @@ -49,6 +49,10 @@ def test_double_excitation_gate(): ] test_list = double_excitation_gate([0, 1, 2, 3]) + for gate in test_list: + if gate.parameters: + print(gate.name, gate.parameters) + # Check gates are correct assert all( control.name == test.name and control.target_qubits == test.target_qubits @@ -69,3 +73,11 @@ def test_givens_excitation_circuit(excitation, expected): control.name == test.name and control.target_qubits == test.target_qubits for control, test in zip(expected, list(test_circuit.queue)) ) + # Check parameters of parametrized gates are all 0.0 + assert all(np.isclose(gate.parameters[0], 0.0) for gate in test_circuit.queue if gate.parameters) + + +def test_givens_excitation_errors(): + """Input excitations are single or double?""" + with pytest.raises(NotImplementedError): + test_circuit = givens_excitation_circuit(4, list(range(6))) From 29fcafc6367db27f5a135ab4d5582692a6cb338f Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 04:21:02 +0000 Subject: [PATCH 18/32] Add convenience function for all excitations --- src/qibochem/ansatz/givens_excitation.py | 70 +++++++++++++++++++++++- 1 file changed, 67 insertions(+), 3 deletions(-) diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py index c0192dc..4a5b6e8 100644 --- a/src/qibochem/ansatz/givens_excitation.py +++ b/src/qibochem/ansatz/givens_excitation.py @@ -4,6 +4,9 @@ from qibo import Circuit, gates +from qibochem.ansatz.hf_reference import hf_circuit +from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations + # Helper functions @@ -81,6 +84,67 @@ def givens_excitation_circuit(n_qubits, excitation, theta=0.0): return circuit -def givens_excitation_ansatz(n_qubits, excitation, n_electrons): - """TODO: Same implementation as UCC?""" - pass +def givens_excitation_ansatz( + molecule, + excitation_level=None, + excitations=None, + include_hf=True, + use_mp2_guess=True, +): + """ + Convenience function for buildng a circuit corresponding to the Givens excitation ansatz with multiple excitations + for a given ``Molecule``. If no excitations are given, it defaults to including all possible spin-allowed + excitations, up to doubles. + + Args: + molecule: The ``Molecule`` of interest. + excitation_level: Include excitations up to how many electrons, i.e. ``"S"`` or ``"D"``. + Ignored if ``excitations`` argument is given. Default: ``"D"``, i.e. double excitations + excitations: List of excitations (e.g. ``[[0, 1, 2, 3], [0, 1, 4, 5]]``) used to build the + UCC circuit. Overrides the ``excitation_level`` argument + include_hf: Whether or not to start the circuit with a Hartree-Fock circuit. Default: ``True`` + use_mp2_guess: Whether to use MP2 amplitudes or a numpy zero array as the initial guess parameter. Default: ``True``; + use the MP2 amplitudes as the default guess parameters + + Returns: + Qibo ``Circuit``: Circuit corresponding to the 'Universal' ansatz + """ + # TODO: Consolidate/Meld this function with the ucc_ansatz function; both are largely identical + + # Get the number of electrons and spin-orbitals from the molecule argument + n_elec = molecule.nelec if molecule.n_active_e is None else molecule.n_active_e + n_orbs = molecule.nso if molecule.n_active_orbs is None else molecule.n_active_orbs + + # Define the excitation level to be used if no excitations given + if excitations is None: + excitation_levels = ("S", "D", "T", "Q") + if excitation_level is None: + excitation_level = "D" + else: + # Check validity of input + assert ( + len(excitation_level) == 1 and excitation_level.upper() in excitation_levels + ), "Unknown input for excitation_level" + # Note: Probably don't be too ambitious and try to do 'T'/'Q' at the moment... + if excitation_level.upper() in ("T", "Q"): + raise NotImplementedError("Cannot handle triple and quadruple excitations!") + # Get the (largest) order of excitation to use + excitation_order = excitation_levels.index(excitation_level.upper()) + 1 + + # Generate and sort all the possible excitations + excitations = [] + for order in range(excitation_order, 0, -1): # Reversed to get higher excitations first + excitations += sort_excitations(generate_excitations(order, range(0, n_elec), range(n_elec, n_orbs))) + else: + # Some checks to ensure the given excitations are valid + assert all(len(_ex) % 2 == 0 for _ex in excitations), "Excitation with an odd number of elements found!" + + # Build the circuit + if include_hf: + circuit = hf_circuit(n_orbs, n_elec) # Only works with (default) JW mapping + else: + circuit = Circuit(n_orbs) + for excitation in excitations: + theta = mp2_amplitude(excitation, molecule.eps, molecule.tei) if use_mp2_guess else None + circuit += givens_excitation_circuit(n_orbs, excitation, theta) + return circuit From b21cc2148a6a2a19ad945963d9f2b745897601c1 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:00:19 +0000 Subject: [PATCH 19/32] Sort out theta argument between two functions `givens_excitation_circuit` and `double_excitation_gate` --- src/qibochem/ansatz/givens_excitation.py | 13 +++++++------ tests/test_givens_excitation.py | 8 ++------ 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py index 4a5b6e8..c44509b 100644 --- a/src/qibochem/ansatz/givens_excitation.py +++ b/src/qibochem/ansatz/givens_excitation.py @@ -10,20 +10,17 @@ # Helper functions -def double_excitation_gate(sorted_orbitals, theta=0.0): +def double_excitation_gate(sorted_orbitals, theta): """ Decomposition of a Givens double excitation gate into single qubit rotations and CNOTs Args: sorted_orbitals (list): Sorted list of orbitals involved in the excitation - theta (float): Rotation angle. Default: 0.0 + theta (float): Rotation angle Returns: (list): List of gates representing the decomposition of the Givens' double excitation gate """ - if theta is None: - theta = 0.0 - result = [] result.append(gates.CNOT(sorted_orbitals[2], sorted_orbitals[3])) result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[2])) @@ -57,7 +54,7 @@ def double_excitation_gate(sorted_orbitals, theta=0.0): # Main function -def givens_excitation_circuit(n_qubits, excitation, theta=0.0): +def givens_excitation_circuit(n_qubits, excitation, theta=None): """ Circuit ansatz for one Givens rotation excitation from Arrazola et al. Reference: https://doi.org/10.22331/q-2022-06-20-742 @@ -66,6 +63,7 @@ def givens_excitation_circuit(n_qubits, excitation, theta=0.0): n_qubits: Number of qubits in the quantum circuit excitation: Iterable of orbitals involved in the excitation; must have an even number of elements E.g. ``[0, 1, 2, 3]`` represents the excitation of electrons in orbitals ``(0, 1)`` to ``(2, 3)`` + theta (float): Rotation angle. Default: 0.0 Returns: Qibo ``Circuit``: Circuit ansatz @@ -74,6 +72,9 @@ def givens_excitation_circuit(n_qubits, excitation, theta=0.0): # Check size of orbitals input assert len(sorted_orbitals) % 2 == 0, f"{excitation} must have an even number of items" + if theta is None: + theta = 0.0 + circuit = Circuit(n_qubits) if len(excitation) == 2: circuit.add(gates.GIVENS(excitation[0], excitation[1], theta)) diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py index 769fa0d..3f7fee1 100644 --- a/tests/test_givens_excitation.py +++ b/tests/test_givens_excitation.py @@ -47,11 +47,7 @@ def test_double_excitation_gate(): gates.CNOT(0, 2), gates.CNOT(2, 3), ] - test_list = double_excitation_gate([0, 1, 2, 3]) - - for gate in test_list: - if gate.parameters: - print(gate.name, gate.parameters) + test_list = double_excitation_gate([0, 1, 2, 3], 0.0) # Check gates are correct assert all( @@ -64,7 +60,7 @@ def test_double_excitation_gate(): "excitation,expected", [ ([0, 2], [gates.GIVENS(0, 2, 0.0)]), - ([0, 1, 2, 3], double_excitation_gate([0, 1, 2, 3])), + ([0, 1, 2, 3], double_excitation_gate([0, 1, 2, 3], 0.0)), ], ) def test_givens_excitation_circuit(excitation, expected): From e55cf97cfbaafdf4c4187108cf3cc757fae7c20c Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:41:13 +0000 Subject: [PATCH 20/32] Simplify ansatz convenience function --- src/qibochem/ansatz/givens_excitation.py | 24 +++--------------------- 1 file changed, 3 insertions(+), 21 deletions(-) diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py index c44509b..5338d9e 100644 --- a/src/qibochem/ansatz/givens_excitation.py +++ b/src/qibochem/ansatz/givens_excitation.py @@ -87,7 +87,6 @@ def givens_excitation_circuit(n_qubits, excitation, theta=None): def givens_excitation_ansatz( molecule, - excitation_level=None, excitations=None, include_hf=True, use_mp2_guess=True, @@ -99,8 +98,6 @@ def givens_excitation_ansatz( Args: molecule: The ``Molecule`` of interest. - excitation_level: Include excitations up to how many electrons, i.e. ``"S"`` or ``"D"``. - Ignored if ``excitations`` argument is given. Default: ``"D"``, i.e. double excitations excitations: List of excitations (e.g. ``[[0, 1, 2, 3], [0, 1, 4, 5]]``) used to build the UCC circuit. Overrides the ``excitation_level`` argument include_hf: Whether or not to start the circuit with a Hartree-Fock circuit. Default: ``True`` @@ -116,29 +113,14 @@ def givens_excitation_ansatz( n_elec = molecule.nelec if molecule.n_active_e is None else molecule.n_active_e n_orbs = molecule.nso if molecule.n_active_orbs is None else molecule.n_active_orbs - # Define the excitation level to be used if no excitations given + # If no excitations given, defaults to all possible double and single excitations if excitations is None: - excitation_levels = ("S", "D", "T", "Q") - if excitation_level is None: - excitation_level = "D" - else: - # Check validity of input - assert ( - len(excitation_level) == 1 and excitation_level.upper() in excitation_levels - ), "Unknown input for excitation_level" - # Note: Probably don't be too ambitious and try to do 'T'/'Q' at the moment... - if excitation_level.upper() in ("T", "Q"): - raise NotImplementedError("Cannot handle triple and quadruple excitations!") - # Get the (largest) order of excitation to use - excitation_order = excitation_levels.index(excitation_level.upper()) + 1 - - # Generate and sort all the possible excitations excitations = [] - for order in range(excitation_order, 0, -1): # Reversed to get higher excitations first + for order in range(2, 0, -1): # Reversed to get double excitations first, then singles excitations += sort_excitations(generate_excitations(order, range(0, n_elec), range(n_elec, n_orbs))) else: # Some checks to ensure the given excitations are valid - assert all(len(_ex) % 2 == 0 for _ex in excitations), "Excitation with an odd number of elements found!" + assert all(len(_ex) in (2, 4) for _ex in excitations), "Only single and double excitations allowed!" # Build the circuit if include_hf: From 17a8cd3715659f36485bf00cbb6c5a153c1dd1f0 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 1 Aug 2024 06:42:26 +0000 Subject: [PATCH 21/32] Add tests for ansatz convenience function TODO: Merge with UCC ansatz convenience function eventually... --- tests/test_givens_excitation.py | 88 +++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py index 3f7fee1..a16ea21 100644 --- a/tests/test_givens_excitation.py +++ b/tests/test_givens_excitation.py @@ -6,11 +6,14 @@ import pytest from qibo import Circuit, gates +from qibochem.ansatz import hf_circuit from qibochem.ansatz.givens_excitation import ( double_excitation_gate, givens_excitation_ansatz, givens_excitation_circuit, ) +from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations +from qibochem.driver import Molecule def test_double_excitation_gate(): @@ -77,3 +80,88 @@ def test_givens_excitation_errors(): """Input excitations are single or double?""" with pytest.raises(NotImplementedError): test_circuit = givens_excitation_circuit(4, list(range(6))) + + +def test_givens_excitation_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))]) + mol.run_pyscf() + + # Build control circuit + control_circuit = hf_circuit(4, 2) + excitations = ([0, 1, 2, 3], [0, 2], [1, 3]) + for excitation in excitations: + control_circuit += givens_excitation_circuit(4, excitation) + + test_circuit = givens_excitation_ansatz(mol) + + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + # for control, test in zip(list(control_circuit.queue), list(test_circuit.queue)) + for control, test in zip(control_circuit.queue, test_circuit.queue) + ) + # Check that number of parametrised gates is the same + assert len(control_circuit.get_parameters()) == len(test_circuit.get_parameters()) + + # Then check that the circuit parameters are the MP2 guess parameters + # Get the MP2 amplitudes first, then expand the list based on the excitation type + mp2_guess_amplitudes = [mp2_amplitude([0, 1, 2, 3], mol.eps, mol.tei) for _ in range(8)] # Doubles + mp2_guess_amplitudes += [0.0, 0.0] # Singles + coeffs = np.array([-0.125, 0.125, -0.125, 0.125, 0.125, -0.125, 0.125, -0.125, 1.0, 1.0]) + mp2_guess_amplitudes = coeffs * np.array(mp2_guess_amplitudes) + # Need to flatten the output of circuit.get_parameters() to compare it to mp2_guess_amplitudes + test_parameters = np.array([_x for _tuple in test_circuit.get_parameters() for _x in _tuple]) + assert np.allclose(mp2_guess_amplitudes, test_parameters) + + +def test_givens_excitation_ansatz_embedding(): + """Test the default arguments of ucc_ansatz using LiH with HF embedding applied, but without the HF circuit""" + mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) + mol.run_pyscf() + mol.hf_embedding(active=[1, 2, 5]) + + # Generate all possible excitations + excitations = [] + for order in range(2, 0, -1): + # 2 electrons, 6 spin-orbitals + excitations += sort_excitations(generate_excitations(order, range(0, 2), range(2, 6))) + # Build control circuit + control_circuit = hf_circuit(6, 0) + for excitation in excitations: + control_circuit += givens_excitation_circuit(6, excitation) + + test_circuit = givens_excitation_ansatz(mol, include_hf=False, use_mp2_guess=False) + + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(list(control_circuit.queue), list(test_circuit.queue)) + ) + # Check that number of parametrised gates is the same + assert len(control_circuit.get_parameters()) == len(test_circuit.get_parameters()) + + # Check that the circuit parameters are all zeros + test_parameters = np.array([_x for _tuple in test_circuit.get_parameters() for _x in _tuple]) + assert np.allclose(test_parameters, np.zeros(len(test_parameters))) + + +def test_ucc_ansatz_excitations(): + """Test the `excitations` argument of ucc_ansatz""" + mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) + mol.run_pyscf() + mol.hf_embedding(active=[1, 2, 5]) + + # Generate all possible excitations + excitations = [[0, 1, 2, 3], [0, 1, 4, 5]] + # Build control circuit + control_circuit = hf_circuit(6, 2) + for excitation in excitations: + control_circuit += givens_excitation_circuit(6, excitation) + + test_circuit = givens_excitation_ansatz(mol, excitations=excitations) + + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(list(control_circuit.queue), list(test_circuit.queue)) + ) + # Check that number of parametrised gates is the same + assert len(control_circuit.get_parameters()) == len(test_circuit.get_parameters()) From ac77c1455823bdfc18aec6748be5dfb22a9e4807 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Sun, 11 Aug 2024 07:46:02 +0000 Subject: [PATCH 22/32] Update and misc. edits to documentation --- README.md | 6 +++- doc/source/api-reference/ansatz.rst | 12 +++++++ doc/source/appendix/citing-qibochem.rst | 40 +++++++++++------------- src/qibochem/ansatz/givens_excitation.py | 9 +++--- src/qibochem/ansatz/symmetry.py | 4 +-- 5 files changed, 42 insertions(+), 29 deletions(-) diff --git a/README.md b/README.md index 69662e9..dc7e398 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,11 @@ If you use the Qibochem plugin please refer to the documentation for citation in ## Contact -For questions, comments and suggestions please contact us at [https://matrix.to/#/#qibo:matrix.org](url) +To get in touch with the community and the developers, consider joining the Qibo workspace on Matrix: + +[![Matrix](https://img.shields.io/matrix/qibo%3Amatrix.org?logo=matrix)](https://matrix.to/#/#qibo:matrix.org) + +If you have a question about the project, contact us at [📫](mailto:qiboteam@qibo.science). ## Contributing diff --git a/doc/source/api-reference/ansatz.rst b/doc/source/api-reference/ansatz.rst index b854440..a271071 100644 --- a/doc/source/api-reference/ansatz.rst +++ b/doc/source/api-reference/ansatz.rst @@ -27,3 +27,15 @@ Basis rotation -------------- .. autofunction:: qibochem.ansatz.basis_rotation.basis_rotation_gates + +Givens Excitation +----------------- + +.. autofunction:: qibochem.ansatz.givens_excitation.givens_excitation_circuit + +.. autofunction:: qibochem.ansatz.givens_excitation.givens_excitation_ansatz + +Symmetry-Preserving +------------------- + +.. autofunction:: qibochem.ansatz.symmetry.symm_preserving_circuit diff --git a/doc/source/appendix/citing-qibochem.rst b/doc/source/appendix/citing-qibochem.rst index b45c97a..1a94e4a 100644 --- a/doc/source/appendix/citing-qibochem.rst +++ b/doc/source/appendix/citing-qibochem.rst @@ -1,22 +1,20 @@ -Publications -============ - -If Qibochem has been significant in your research, and you would like to acknowledge -the project in your academic publication, we suggest citing the following documents: - -Software References in Zenodo ------------------------------ - -Wong Zicheng, Adrian Mak, Tan Le, Stefano Carrazza, Alessandro Candido, & Ye Jun. (2024). qiboteam/qibochem: Qibochem 0.0.1 (v0.0.1). Zenodo. https://doi.org/10.5281/zenodo.10473173 - -Authorship Guideline --------------------- - -In order to appear as an author of a Qibochem publication (paper, proceedings, etc) -each author must fullfil the following requirements: - -* Participate to the official meetings. - -* Contribute to the code with documented commits. - +Publications +============ + +If Qibochem has been significant in your research, and you would like to acknowledge the project in your academic publication, we suggest citing the following reference: + +Software References in Zenodo +----------------------------- + +Wong Zi Cheng, Adrian Matthew Mak, Tan Le, Stefano Carrazza, Alessandro Candido & Ye Jun. (2024). qiboteam/qibochem: Qibochem 0.0.1 (v0.0.1). Zenodo. https://doi.org/10.5281/zenodo.10473173 + +Authorship Guideline +-------------------- + +In order to appear as an author of a Qibochem publication (paper, proceedings, etc), one must fullfil the following requirements: + +* Participate to the official meetings. + +* Contribute to the code with documented commits. + * Contribute to the manuscript elaboration. diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py index 5338d9e..3d8c3ef 100644 --- a/src/qibochem/ansatz/givens_excitation.py +++ b/src/qibochem/ansatz/givens_excitation.py @@ -56,17 +56,16 @@ def double_excitation_gate(sorted_orbitals, theta): # Main function def givens_excitation_circuit(n_qubits, excitation, theta=None): """ - Circuit ansatz for one Givens rotation excitation from Arrazola et al. Reference: - https://doi.org/10.22331/q-2022-06-20-742 + Circuit ansatz corresponding to the Givens rotation excitation from Arrazola et al. (https://doi.org/10.22331/q-2022-06-20-742) for a single excitation. Args: - n_qubits: Number of qubits in the quantum circuit + n_qubits: Number of qubits in the circuit excitation: Iterable of orbitals involved in the excitation; must have an even number of elements E.g. ``[0, 1, 2, 3]`` represents the excitation of electrons in orbitals ``(0, 1)`` to ``(2, 3)`` theta (float): Rotation angle. Default: 0.0 Returns: - Qibo ``Circuit``: Circuit ansatz + Qibo ``Circuit``: Circuit ansatz for a single Givens excitation """ sorted_orbitals = sorted(excitation) # Check size of orbitals input @@ -105,7 +104,7 @@ def givens_excitation_ansatz( use the MP2 amplitudes as the default guess parameters Returns: - Qibo ``Circuit``: Circuit corresponding to the 'Universal' ansatz + Qibo ``Circuit``: Circuit corresponding to a Givens excitation circuit ansatz """ # TODO: Consolidate/Meld this function with the ucc_ansatz function; both are largely identical diff --git a/src/qibochem/ansatz/symmetry.py b/src/qibochem/ansatz/symmetry.py index 170df7b..2a4992e 100644 --- a/src/qibochem/ansatz/symmetry.py +++ b/src/qibochem/ansatz/symmetry.py @@ -85,14 +85,14 @@ def a_gate_indices(n_qubits, n_electrons, x_gates): # Main function def symm_preserving_circuit(n_qubits, n_electrons): """ - Symmetry-preserving circuit ansatz from Gard et al. Reference: https://doi.org/10.1038/s41534-019-0240-1 + Symmetry-preserving circuit ansatz from Gard et al. (https://doi.org/10.1038/s41534-019-0240-1) Args: n_qubits: Number of qubits in the quantum circuit n_electrons: Number of electrons in the molecular system Returns: - Qibo ``Circuit``: Circuit ansatz + Qibo ``Circuit``: Circuit corresponding to the symmetry-preserving ansatz """ circuit = Circuit(n_qubits) x_gates = x_gate_indices(n_qubits, n_electrons) From f1cf698effb364438b1670277382d4a72f65182f Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Sat, 17 Aug 2024 06:15:06 +0000 Subject: [PATCH 23/32] Start writing adaptive example in tutorial --- doc/source/tutorials/index.rst | 1 + src/qibochem/ansatz/__init__.py | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/doc/source/tutorials/index.rst b/doc/source/tutorials/index.rst index 4ce96f7..3617a44 100644 --- a/doc/source/tutorials/index.rst +++ b/doc/source/tutorials/index.rst @@ -12,3 +12,4 @@ This section provides examples of how to use Qibochem. hamiltonian ansatz measurement + adaptive diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py index e0b77f5..00609da 100644 --- a/src/qibochem/ansatz/__init__.py +++ b/src/qibochem/ansatz/__init__.py @@ -1,4 +1,8 @@ from qibochem.ansatz.basis_rotation import basis_rotation_gates +from qibochem.ansatz.givens_excitation import ( + givens_excitation_ansatz, + givens_excitation_circuit, +) from qibochem.ansatz.hardware_efficient import he_circuit from qibochem.ansatz.hf_reference import hf_circuit from qibochem.ansatz.symmetry import symm_preserving_circuit From 5a9acd4e2e836cca6e1b4e3183b521eb3396beb7 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Tue, 20 Aug 2024 07:21:00 +0000 Subject: [PATCH 24/32] Update code for single excitation --- src/qibochem/ansatz/givens_excitation.py | 29 +++++++++++++++++++++--- tests/test_givens_excitation.py | 28 ++++++++++++++++++++--- 2 files changed, 51 insertions(+), 6 deletions(-) diff --git a/src/qibochem/ansatz/givens_excitation.py b/src/qibochem/ansatz/givens_excitation.py index 3d8c3ef..e891d15 100644 --- a/src/qibochem/ansatz/givens_excitation.py +++ b/src/qibochem/ansatz/givens_excitation.py @@ -1,5 +1,6 @@ """ -Module documentation +Circuit ansatz for representing a fermionic excitation as a Givens rotation by Arrazola et al. +Reference: https://doi.org/10.22331/q-2022-06-20-742 """ from qibo import Circuit, gates @@ -10,6 +11,28 @@ # Helper functions +def single_excitation_gate(sorted_orbitals, theta): + """ + Decomposition of a Givens single excitation gate into single qubit rotations and CNOTs. In principle, should be + identical to gates.GIVENS(qubit0, qubit1, theta) + + Args: + sorted_orbitals (list): Sorted list of orbitals involved in the excitation + theta (float): Rotation angle + + Returns: + (list): List of gates representing the decomposition of the Givens' single excitation gate + """ + result = [] + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1])) + result.append(gates.RY(sorted_orbitals[0], 0.5 * theta)) + result.append(gates.CNOT(sorted_orbitals[1], sorted_orbitals[0])) + result.append(gates.RY(sorted_orbitals[0], -0.5 * theta)) + result.append(gates.CNOT(sorted_orbitals[1], sorted_orbitals[0])) + result.append(gates.CNOT(sorted_orbitals[0], sorted_orbitals[1])) + return result + + def double_excitation_gate(sorted_orbitals, theta): """ Decomposition of a Givens double excitation gate into single qubit rotations and CNOTs @@ -56,7 +79,7 @@ def double_excitation_gate(sorted_orbitals, theta): # Main function def givens_excitation_circuit(n_qubits, excitation, theta=None): """ - Circuit ansatz corresponding to the Givens rotation excitation from Arrazola et al. (https://doi.org/10.22331/q-2022-06-20-742) for a single excitation. + Circuit ansatz corresponding to the Givens rotation excitation from Arrazola et al. Args: n_qubits: Number of qubits in the circuit @@ -76,7 +99,7 @@ def givens_excitation_circuit(n_qubits, excitation, theta=None): circuit = Circuit(n_qubits) if len(excitation) == 2: - circuit.add(gates.GIVENS(excitation[0], excitation[1], theta)) + circuit.add(single_excitation_gate(sorted_orbitals, theta)) elif len(excitation) == 4: circuit.add(double_excitation_gate(sorted_orbitals, theta)) else: diff --git a/tests/test_givens_excitation.py b/tests/test_givens_excitation.py index a16ea21..2344a20 100644 --- a/tests/test_givens_excitation.py +++ b/tests/test_givens_excitation.py @@ -11,11 +11,33 @@ double_excitation_gate, givens_excitation_ansatz, givens_excitation_circuit, + single_excitation_gate, ) from qibochem.ansatz.util import generate_excitations, mp2_amplitude, sort_excitations from qibochem.driver import Molecule +def test_single_excitation_gate(): + # Hardcoded test + theta = 0.1 + + control_gates = [ + gates.CNOT(0, 1), + gates.RY(0, 0.5 * theta), + gates.CNOT(1, 0), + gates.RY(0, -0.5 * theta), + gates.CNOT(1, 0), + gates.CNOT(0, 1), + ] + test_list = single_excitation_gate([0, 1, 2, 3], 0.1) + + # Check gates are correct + assert all( + control.name == test.name and control.target_qubits == test.target_qubits + for control, test in zip(control_gates, test_list) + ) + + def test_double_excitation_gate(): # Hardcoded test theta = 0.0 @@ -62,7 +84,7 @@ def test_double_excitation_gate(): @pytest.mark.parametrize( "excitation,expected", [ - ([0, 2], [gates.GIVENS(0, 2, 0.0)]), + ([0, 2], single_excitation_gate([0, 2], 0.0)), ([0, 1, 2, 3], double_excitation_gate([0, 1, 2, 3], 0.0)), ], ) @@ -106,8 +128,8 @@ def test_givens_excitation_ansatz_h2(): # Then check that the circuit parameters are the MP2 guess parameters # Get the MP2 amplitudes first, then expand the list based on the excitation type mp2_guess_amplitudes = [mp2_amplitude([0, 1, 2, 3], mol.eps, mol.tei) for _ in range(8)] # Doubles - mp2_guess_amplitudes += [0.0, 0.0] # Singles - coeffs = np.array([-0.125, 0.125, -0.125, 0.125, 0.125, -0.125, 0.125, -0.125, 1.0, 1.0]) + mp2_guess_amplitudes += [0.0, 0.0, 0.0, 0.0] # Singles + coeffs = np.array([-0.125, 0.125, -0.125, 0.125, 0.125, -0.125, 0.125, -0.125, 1.0, 1.0, 1.0, 1.0]) mp2_guess_amplitudes = coeffs * np.array(mp2_guess_amplitudes) # Need to flatten the output of circuit.get_parameters() to compare it to mp2_guess_amplitudes test_parameters = np.array([_x for _tuple in test_circuit.get_parameters() for _x in _tuple]) From 120b63f13a0d30503351b4eb7c71628402d82287 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 21 Aug 2024 01:57:14 +0000 Subject: [PATCH 25/32] Draft tutorial for adaptive methods --- doc/source/tutorials/adaptive.rst | 357 ++++++++++++++++++++++++++++++ 1 file changed, 357 insertions(+) create mode 100644 doc/source/tutorials/adaptive.rst diff --git a/doc/source/tutorials/adaptive.rst b/doc/source/tutorials/adaptive.rst new file mode 100644 index 0000000..ea58bce --- /dev/null +++ b/doc/source/tutorials/adaptive.rst @@ -0,0 +1,357 @@ +Adaptive methods +================ + +Despite VQE being cheaper than QPE, circuit depth is still a big problem for today's quantum hardware. + +.. code-block:: python + + from qibochem.driver import Molecule + from qibochem.ansatz import ucc_ansatz + + mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) + mol.run_pyscf() + circuit = ucc_ansatz(mol) + print(circuit.summary()) + +Output: + +.. code-block:: output + + Circuit depth = 8937 + Total number of gates = 15108 + Number of qubits = 12 + Most common gates: + cx: 6976 + h: 4992 + sdg: 1248 + s: 1248 + rz: 640 + x: 4 + + +As shown in the above code block, a full UCCSD circuit for a LiH/STO-3G system has a circuit depth of 8937 (!), with almost 7000 CNOT gates required! +Hence, there is still a need to further reduce and simplify the circuit ansatzes used for running a VQE simulation. + +Other than designing shorter circuit ansatzes, one alternative approach is through the use of energy gradients - for instance, through the Parameter Shift Rule on hardware - to filter and reduce the number of fermionic excitations in a circuit ansatz. (REFS!!!) +This is known as an adaptive method, in the sense that the quantum gates used to construct the circuit ansatz, as well as its actual structure and arrangement is not fixed, and varies depending on the molecular system under study. + +For example, in a H2/STO-3G system mapped with the Jordan-Wigner transformation, there are three possible spin-allowed fermionic excitations: +two single excitations (``[0, 2]``, ``[1, 3]``) and one double excitation (``[0, 1, 2, 3]``). +The full UCCSD circuit for this system has been shown in an earlier example (ref), and it requires 64 CNOT gates for this simple molecular system. + +Let's look at the gradients of each of the individual fermionic excitations: + +.. code-block:: python + + from qibo.derivative import parameter_shift + + from qibochem.driver import Molecule + from qibochem.ansatz import hf_circuit, ucc_circuit + + mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) + mol.run_pyscf() + hamiltonian = mol.hamiltonian() + + excitations = [[0, 1, 2, 3], [0, 2], [1, 3]] + circuit = hf_circuit(4, 2) + for excitation in excitations: + _circuit = circuit.copy() + _circuit += ucc_circuit(4, excitation) + n_parameters = len(_circuit.get_parameters()) + gradients = [round(parameter_shift(_circuit, hamiltonian, index), 5) for index in range(n_parameters)] + print(f"Energy gradients for {excitation}: {gradients}") + +Output: + +.. code-block:: output + + Energy gradients for [0, 1, 2, 3]: [0.179, -0.179, -0.179, -0.179, 0.179, 0.179, 0.179, -0.179] + Energy gradients for [0, 2]: [0.0, 0.0] + Energy gradients for [1, 3]: [0.0, 0.0] + +Only the doubles excitation has a non-zero energy gradient, which follows Brillouin's theorem. +Running the VQE with only the doubles excitation then gives: + +.. code-block:: python + + import numpy as np + + from qibo.models import VQE + + from qibochem.driver import Molecule + from qibochem.ansatz import hf_circuit, ucc_circuit + + mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) + mol.run_pyscf() + hamiltonian = mol.hamiltonian() + + circuit = hf_circuit(4, 2) + circuit += ucc_circuit(4, [0, 1, 2, 3]) + + vqe = VQE(circuit, hamiltonian) + + # Optimize starting from a random guess for the variational parameters + initial_parameters = np.random.uniform(0, 2*np.pi, len(circuit.get_parameters())) + best, params, extra = vqe.minimize(initial_parameters, method='BFGS', compile=False) + + # Exact result + print(f"Exact result: {mol.eigenvalues(hamiltonian)[0]:.7f}") + print(f" VQE result: {best:.7f}") + +Output: + +.. code-block:: output + + Exact result: -1.1361895 + VQE result: -1.1361895 + +As can be seen, we managed to find the exact result using only one doubles excitation! + +Next, let's look at the potential savings for the LiH/STO-3G system. +To reduce the circuit depth further, we will use the more modern ansatz, the Givens excitation circuit from Arrazola et al., instead of the UCC ansatz. + + + +TODO: +Starting with a HF ansatz, find gradient of each Givens excitation, and sort the excitations by the magnitude of the gradient. + +Then apply circuit ansatz for each excitation iteratively, until energy converges. +How much time/circuit depth do we save in this approach, compared to the naive, add everything approach? + + +.. comment + OLD STUFF + --------- + + A quantum circuit comprising parameterized gates (`e.g.` :math:`RX(\theta)`, :math:`RY(\theta)` and :math:`RZ(\theta)`), + represents a unitary transformation :math:`U(\theta)` that transforms some initial quantum state into a parametrized ansatz state :math:`|\psi(\theta)\rangle`. + + Examples of some ansatzes available in Qibochem are described in the subsections below. + + Hardware Efficient Ansatz + ------------------------- + + Qibochem provides a hardware efficient ansatz that simply consists of a layer of single-qubit rotation gates followed by a layer of two-qubit gates that entangle the qubits. + For the H\ :sub:`2` case discussed in previous sections, a possible hardware efficient circuit ansatz can be constructed as such: + + .. image:: qibochem_doc_ansatz_hardware-efficient.svg + + .. code-block:: python + + from qibochem.ansatz import he_circuit + + nqubits = 4 + nlayers = 1 + + circuit = he_circuit(nqubits, nlayers) + print(circuit.draw()) + + .. code-block:: output + + q0: ─RY─RZ─o─────Z─ + q1: ─RY─RZ─Z─o───|─ + q2: ─RY─RZ───Z─o─|─ + q3: ─RY─RZ─────Z─o─ + + The energy of the state generated from the hardware efficient ansatz for the fermionic two-body Hamiltonian can then be estimated, using state vectors or samples. + + The following example demonstrates how the energy of the H2 molecule is affected with respect to the rotational parameters: + + .. code-block:: python + + import numpy as np + + from qibochem.driver import Molecule + from qibochem.measurement.expectation import expectation + 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() + hamiltonian = mol.hamiltonian() + + # Define and build the HEA + nlayers = 1 + nqubits = mol.nso + ntheta = 2 * nqubits * nlayers + circuit = he_circuit(nqubits, nlayers) + + print("Energy expectation values for thetas: ") + print("-----------------------------") + print("| theta | Electronic energy |") + print("|---------------------------|") + thetas = [-0.2, 0.0, 0.2] + for theta in thetas: + params = np.full(ntheta, theta) + circuit.set_parameters(params) + electronic_energy = expectation(circuit, hamiltonian) + print(f"| {theta:5.1f} | {electronic_energy:^18.12f}|") + print("-----------------------------") + + + .. code-block:: output + + converged SCF energy = -1.11628373627429 + + Energy expectation values for thetas: + ----------------------------- + | theta | Electronic energy | + |---------------------------| + | -0.2 | 0.673325849299 | + | 0.0 | 0.707418334474 | + | 0.2 | 0.673325849299 | + ----------------------------- + + + .. _UCC Ansatz: + + Unitary Coupled Cluster Ansatz + ------------------------------ + + The Unitary Coupled Cluster (UCC) ansatz [#f1]_ [#f2]_ [#f3]_ is a variant of the popular gold standard Coupled Cluster ansatz [#f4]_ of quantum chemistry. + The UCC wave function is a parameterized unitary transformation of a reference wave function :math:`\psi_{\mathrm{ref}}`, of which a common choice is the Hartree-Fock wave function. + + .. math:: + + \begin{align*} + |\psi_{\mathrm{UCC}}\rangle &= U(\theta)|\psi_{\mathrm{ref}}\rangle \\ + &= e^{\hat{T}(\theta) - \hat{T}^\dagger(\theta)}|\psi_{\mathrm{ref}}\rangle + \end{align*} + + + Similar to the process for the molecular Hamiltonian, the fermionic excitation operators :math:`\hat{T}` and :math:`\hat{T}^\dagger` are mapped using e.g. Jordan-Wigner mapping into Pauli operators. + This is typically followed by a Suzuki-Trotter decomposition of the exponentials of these Pauli operators, which allows the UCC ansatz to be implemented on quantum computers. [#f5]_ + + An example of how to build a UCC doubles circuit ansatz for the :math:`H_2` molecule is given as: + + .. code-block:: python + + from qibochem.driver.molecule import Molecule + 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() + hamiltonian = mol.hamiltonian() + + # Set parameters for the rest of the experiment + n_qubits = mol.nso + n_electrons = mol.nelec + + # Build UCCD circuit + circuit = hf_circuit(n_qubits, n_electrons) # Start with HF circuit + circuit += ucc_circuit(n_qubits, [0, 1, 2, 3]) # Then add the double excitation circuit ansatz + + print(circuit.draw()) + + .. code-block:: output + + q0: ─X──H─────X─RZ─X─────H──RX─────X─RZ─X─────RX─RX─────X─RZ─X─────RX─H─── ... + q1: ─X──H───X─o────o─X───H──RX───X─o────o─X───RX─H────X─o────o─X───H──RX── ... + q2: ─RX───X─o────────o─X─RX─RX─X─o────────o─X─RX─H──X─o────────o─X─H──H──X ... + q3: ─H────o────────────o─H──H──o────────────o─H──H──o────────────o─H──H──o ... + + q0: ... ───X─RZ─X─────H──RX─────X─RZ─X─────RX─H──────X─RZ─X─────H──H──────X─RZ ... + q1: ... ─X─o────o─X───RX─H────X─o────o─X───H──RX───X─o────o─X───RX─H────X─o─── ... + q2: ... ─o────────o─X─H──RX─X─o────────o─X─RX─RX─X─o────────o─X─RX─H──X─o───── ... + q3: ... ────────────o─H──RX─o────────────o─RX─RX─o────────────o─RX─RX─o─────── ... + + q0: ... ─X─────H──RX─────X─RZ─X─────RX─ + q1: ... ─o─X───H──RX───X─o────o─X───RX─ + q2: ... ───o─X─H──H──X─o────────o─X─H── + q3: ... ─────o─RX─RX─o────────────o─RX─ + + + .. + _Basis rotation ansatz + + Basis rotation ansatz + --------------------- + + The starting points for contemporary quantum chemistry methods are often those based on the mean field approximation within a (finite) molecular orbital basis, i.e. the Hartree-Fock method. The electronic energy is calculated as the mean value of the electronic Hamiltonian :math:`\hat{H}_{\mathrm{elec}}` acting on a normalized single Slater determinant function :math:`\psi` [#f6]_ + + .. math:: + + \begin{align*} + E[\psi] &= \langle \psi | \hat{H}_{\mathrm{elec}} |\psi \rangle \\ + &= \sum_i^{N_f} \langle \phi_i |\hat{h}|\phi_i \rangle + \frac{1}{2} \sum_{i,j}^{N_f} + \langle \phi_i\phi_j||\phi_i\phi_j \rangle + \end{align*} + + The orthonormal molecular orbitals :math:`\phi` are optimized by a direct minimization of the energy functional with respect to parameters :math:`\kappa` that parameterize the unitary rotations of the orbital basis. Qibochem's implementation uses the QR decomposition of the unitary matrix as employed by Clements et al., [#f7]_ which results in a rectangular gate layout of `Givens rotation gates `_ that yield linear CNOT gate depth when decomposed. + + + .. code-block:: python + + import numpy as np + from qibochem.driver.molecule import Molecule + from qibochem.ansatz import basis_rotation, ucc + from qibo import Circuit, gates, models + + def basis_rotation_circuit(mol, parameters=0.0): + + nqubits = mol.nso + occ = range(0, mol.nelec) + vir = range(mol.nelec, mol.nso) + + U, kappa = basis_rotation.unitary(occ, vir, parameters=parameters) + gate_angles, final_U = basis_rotation.givens_qr_decompose(U) + gate_layout = basis_rotation.basis_rotation_layout(nqubits) + gate_list, ordered_angles = basis_rotation.basis_rotation_gates(gate_layout, gate_angles, kappa) + + circuit = Circuit(nqubits) + for _i in range(mol.nelec): + circuit.add(gates.X(_i)) + circuit.add(gate_list) + + return circuit, gate_angles + + h3p = Molecule([('H', (0.0000, 0.0000, 0.0000)), + ('H', (0.0000, 0.0000, 0.8000)), + ('H', (0.0000, 0.0000, 1.6000))], + charge=1, multiplicity=1) + h3p.run_pyscf(max_scf_cycles=1) + + e_init = h3p.e_hf + h3p_sym_ham = h3p.hamiltonian("sym", h3p.oei, h3p.tei, 0.0, "jw") + + hf_circuit, qubit_parameters = basis_rotation_circuit(h3p, parameters=0.1) + + print(hf_circuit.draw()) + + vqe = models.VQE(hf_circuit, h3p_sym_ham) + res = vqe.minimize(qubit_parameters) + + print('energy of initial guess: ', e_init) + print('energy after vqe : ', res[0]) + + .. code-block:: output + + q0: ─X─G─────────G─────────G───────── + q1: ─X─G─────G───G─────G───G─────G─── + q2: ─────G───G─────G───G─────G───G─── + q3: ─────G─────G───G─────G───G─────G─ + q4: ───────G───G─────G───G─────G───G─ + q5: ───────G─────────G─────────G───── + basis rotation: using uniform value of 0.1 for each parameter value + energy of initial guess: -1.1977713400022736 + energy after vqe : -1.2024564133305427 + + + + + + + .. rubric:: References + + .. [#f1] Kutzelnigg, W. (1977). 'Pair Correlation Theories', in Schaefer, H.F. (eds) Methods of Electronic Structure Theory. Modern Theoretical Chemistry, vol 3. Springer, Boston, MA. + + .. [#f2] Whitfield, J. D. et al., 'Simulation of Electronic Structure Hamiltonians using Quantum Computers', Mol. Phys. 109 (2011) 735. + + .. [#f3] Anand. A. et al., 'A Quantum Computing view on Unitary Coupled Cluster Theory', Chem. Soc. Rev. 51 (2022) 1659. + + .. [#f4] Crawford, T. D. et al., 'An Introduction to Coupled Cluster Theory for Computational Chemists', in Reviews in Computational Chemistry 14 (2007) 33. + + .. [#f5] Barkoutsos, P. K. et al., 'Quantum algorithms for electronic structure calculations: Particle-hole Hamiltonian and optimized wave-function expansions', Phys. Rev. A 98 (2018) 022322. + + .. [#f6] Piela, L. (2007). 'Ideas of Quantum Chemistry'. Elsevier B. V., the Netherlands. + + .. [#f7] Clements, W. R. et al., 'Optimal Design for Universal Multiport Interferometers', Optica 3 (2016) 1460. From 32f65ce7f208f259f148dbe3bb9c92551cf7f072 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Wed, 21 Aug 2024 10:07:33 +0000 Subject: [PATCH 26/32] Add more to adaptive tutorial --- doc/source/tutorials/adaptive.rst | 111 +++++++++++++++++++++++++----- 1 file changed, 95 insertions(+), 16 deletions(-) diff --git a/doc/source/tutorials/adaptive.rst b/doc/source/tutorials/adaptive.rst index ea58bce..27538de 100644 --- a/doc/source/tutorials/adaptive.rst +++ b/doc/source/tutorials/adaptive.rst @@ -10,6 +10,7 @@ Despite VQE being cheaper than QPE, circuit depth is still a big problem for tod mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) mol.run_pyscf() + mol.hf_embedding(active=[0, 1, 2, 5]) circuit = ucc_ansatz(mol) print(circuit.summary()) @@ -17,22 +18,21 @@ Output: .. code-block:: output - Circuit depth = 8937 - Total number of gates = 15108 - Number of qubits = 12 + Circuit depth = 1874 + Total number of gates = 3300 + Number of qubits = 8 Most common gates: - cx: 6976 - h: 4992 - sdg: 1248 - s: 1248 - rz: 640 + cx: 1312 + h: 1216 + sdg: 304 + s: 304 + rz: 160 x: 4 - -As shown in the above code block, a full UCCSD circuit for a LiH/STO-3G system has a circuit depth of 8937 (!), with almost 7000 CNOT gates required! +As shown in the above code block, the full UCCSD circuit for a simplified LiH/STO-3G system has a circuit depth of 1874 (!), with more than 1000 CNOT gates required! Hence, there is still a need to further reduce and simplify the circuit ansatzes used for running a VQE simulation. -Other than designing shorter circuit ansatzes, one alternative approach is through the use of energy gradients - for instance, through the Parameter Shift Rule on hardware - to filter and reduce the number of fermionic excitations in a circuit ansatz. (REFS!!!) +Other than designing shorter and more efficient circuit ansatzes, one alternative approach is through the use of energy gradients - for instance, through the Parameter Shift Rule on hardware - to filter and reduce the number of fermionic excitations in a circuit ansatz. (REFS!!!) This is known as an adaptive method, in the sense that the quantum gates used to construct the circuit ansatz, as well as its actual structure and arrangement is not fixed, and varies depending on the molecular system under study. For example, in a H2/STO-3G system mapped with the Jordan-Wigner transformation, there are three possible spin-allowed fermionic excitations: @@ -105,17 +105,96 @@ Output: Exact result: -1.1361895 VQE result: -1.1361895 -As can be seen, we managed to find the exact result using only one doubles excitation! +We managed to find the exact result by applying only the doubles excitation! -Next, let's look at the potential savings for the LiH/STO-3G system. +Next, let's look at the potential savings for the simplified LiH/STO-3G system. To reduce the circuit depth further, we will use the more modern ansatz, the Givens excitation circuit from Arrazola et al., instead of the UCC ansatz. +As was done in the above example, we will start with a HF circuit, then find the gradients for each circuit ansatz corresponding to a fermionic excitation. +After that, the excitation with the largest absolute value of the gradient will be added to the initial circuit, followed by a VQE simulation. + +.. code-block:: python + + from qibo.derivative import parameter_shift + from qibo.models import VQE + + from qibochem.driver import Molecule + from qibochem.ansatz import hf_circuit, givens_excitation_circuit, generate_excitations, sort_excitations + + mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) + mol.run_pyscf() + mol.hf_embedding(active=[0, 1, 2, 5]) + hamiltonian = mol.hamiltonian() + + n_qubits = mol.n_active_orbs + n_elec = mol.n_active_e + circuit = hf_circuit(n_qubits, n_elec) -TODO: -Starting with a HF ansatz, find gradient of each Givens excitation, and sort the excitations by the magnitude of the gradient. + excitations = sort_excitations(generate_excitations(2, list(range(n_elec)), list(range(n_elec, n_qubits)))) + excitations += sort_excitations(generate_excitations(1, list(range(n_elec)), list(range(n_elec, n_qubits)))) + + excitation_gradients = {} + for excitation in excitations: + _circuit = circuit.copy() + _circuit += givens_excitation_circuit(n_qubits, excitation) + n_parameters = len(_circuit.get_parameters()) + gradient = [round(parameter_shift(_circuit, hamiltonian, index), 5) for index in range(n_parameters)] + print(f"Energy gradients for {excitation}: {gradient}") + excitation_gradients[tuple(excitation)] = gradient[0] # Gradient magnitude is equal throughout + + # Find the excitation corresponding to the largest gradient, and add it to the circuit + max_grad = max(excitation_gradients, key=lambda x: abs(excitation_gradientis.get(x))) + print(f"\nExcitation with the largest gradient: {max_grad}; Gradient = {excitation_gradients[max_grad]}") + circuit += givens_excitation_circuit(n_qubits, max_grad) + + # Run VQE with the updated circuit + vqe = VQE(circuit, hamiltonian) + + circuit_parameters = [param for _tuple in circuit.get_parameters() for param in _tuple] + best, params, extra = vqe.minimize(circuit_parameters, method='BFGS', compile=False) + + print(f" HF energy: {mol.e_hf:.7f}") + print(f"VQE result: {best:.7f}") + +Output: + +.. code-block:: output -Then apply circuit ansatz for each excitation iteratively, until energy converges. + Energy gradients for [0, 1, 4, 5]: [0.02132, -0.02132, 0.02132, -0.02132, -0.02132, 0.02132, -0.02132, 0.02132] + Energy gradients for [0, 1, 6, 7]: [0.00569, -0.00569, 0.00569, -0.00569, -0.00569, 0.00569, -0.00569, 0.00569] + Energy gradients for [2, 3, 4, 5]: [0.01136, -0.01136, 0.01136, -0.01136, -0.01136, 0.01136, -0.01136, 0.01136] + Energy gradients for [2, 3, 6, 7]: [0.12225, -0.12225, 0.12225, -0.12225, -0.12225, 0.12225, -0.12225, 0.12225] + Energy gradients for [0, 1, 4, 7]: [0.00016, -0.00016, 0.00016, -0.00016, -0.00016, 0.00016, -0.00016, 0.00016] + Energy gradients for [0, 1, 5, 6]: [-0.00016, 0.00016, -0.00016, 0.00016, 0.00016, -0.00016, 0.00016, -0.00016] + Energy gradients for [2, 3, 4, 7]: [-0.03254, 0.03254, -0.03254, 0.03254, 0.03254, -0.03254, 0.03254, -0.03254] + Energy gradients for [2, 3, 5, 6]: [0.03254, -0.03254, 0.03254, -0.03254, -0.03254, 0.03254, -0.03254, 0.03254] + Energy gradients for [0, 3, 4, 5]: [0.00029, -0.00029, 0.00029, -0.00029, -0.00029, 0.00029, -0.00029, 0.00029] + Energy gradients for [1, 2, 4, 5]: [-0.00029, 0.00029, -0.00029, 0.00029, 0.00029, -0.00029, 0.00029, -0.00029] + Energy gradients for [0, 3, 6, 7]: [0.00108, -0.00108, 0.00108, -0.00108, -0.00108, 0.00108, -0.00108, 0.00108] + Energy gradients for [1, 2, 6, 7]: [-0.00108, 0.00108, -0.00108, 0.00108, 0.00108, -0.00108, 0.00108, -0.00108] + Energy gradients for [0, 2, 4, 6]: [0.00299, -0.00299, 0.00299, -0.00299, -0.00299, 0.00299, -0.00299, 0.00299] + Energy gradients for [1, 3, 5, 7]: [0.00299, -0.00299, 0.00299, -0.00299, -0.00299, 0.00299, -0.00299, 0.00299] + Energy gradients for [0, 3, 4, 7]: [-0.00236, 0.00236, -0.00236, 0.00236, 0.00236, -0.00236, 0.00236, -0.00236] + Energy gradients for [0, 3, 5, 6]: [-0.00063, 0.00063, -0.00063, 0.00063, 0.00063, -0.00063, 0.00063, -0.00063] + Energy gradients for [1, 2, 4, 7]: [-0.00063, 0.00063, -0.00063, 0.00063, 0.00063, -0.00063, 0.00063, -0.00063] + Energy gradients for [1, 2, 5, 6]: [-0.00236, 0.00236, -0.00236, 0.00236, 0.00236, -0.00236, 0.00236, -0.00236] + Energy gradients for [0, 4]: [0.0, -0.0] + Energy gradients for [1, 5]: [-0.0, 0.0] + Energy gradients for [0, 6]: [0.0, -0.0] + Energy gradients for [1, 7]: [-0.0, 0.0] + Energy gradients for [2, 4]: [-0.0, 0.0] + Energy gradients for [3, 5]: [0.0, -0.0] + Energy gradients for [2, 6]: [-0.0, 0.0] + Energy gradients for [3, 7]: [0.0, -0.0] + + Excitation with the largest gradient: (2, 3, 6, 7); Gradient = 0.12225 + HF energy: -7.8605387 + VQE result: -7.8732886 + + +Energy difference of ~0.01 Hartrees, result still yet to converge. +So apply circuit ansatz for each excitation with the largest gradient iteratively, until energy converges. How much time/circuit depth do we save in this approach, compared to the naive, add everything approach? From a928d1f8417e7045badccdc3701ce92ad1534e59 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 22 Aug 2024 08:46:24 +0000 Subject: [PATCH 27/32] Finish of text of adaptive tutorial Left with cleaning up and references --- doc/source/tutorials/adaptive.rst | 167 ++++++++++++++++++++++++++- doc/source/tutorials/measurement.rst | 5 +- 2 files changed, 168 insertions(+), 4 deletions(-) diff --git a/doc/source/tutorials/adaptive.rst b/doc/source/tutorials/adaptive.rst index 27538de..9dbf19e 100644 --- a/doc/source/tutorials/adaptive.rst +++ b/doc/source/tutorials/adaptive.rst @@ -193,9 +193,170 @@ Output: VQE result: -7.8732886 -Energy difference of ~0.01 Hartrees, result still yet to converge. -So apply circuit ansatz for each excitation with the largest gradient iteratively, until energy converges. -How much time/circuit depth do we save in this approach, compared to the naive, add everything approach? +After adding the circuit ansatz corresponding to one double excitation and running VQE, +the resultant energy was found to be about 0.01 Hartrees lower compared to the bare HF ansatz. +Clearly, there is still room for further improvement in the obtained energies. + +To do so, we further extend the circuit ansatz by adding in more excitations iteratively, +with the excitation with the largest (in magnitute) energy gradients added on at each step. +This can be carried out until the difference between each iteration is small (<0.001 Hartrees), or until there are no more remaining excitations to be added. + +.. warning:: + + The code block below might take a few minutes to run! + +.. code-block:: python + + import numpy as np + + from qibo.models import VQE + from qibo.derivative import parameter_shift + + from qibochem.driver import Molecule + from qibochem.ansatz import hf_circuit, givens_excitation_circuit, generate_excitations, sort_excitations + + mol = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) + mol.run_pyscf() + mol.hf_embedding(active=[0, 1, 2, 5]) + hamiltonian = mol.hamiltonian() + + exact_result = mol.eigenvalues(hamiltonian)[0] + + n_qubits = mol.n_active_orbs + n_elec = mol.n_active_e + + circuit = hf_circuit(n_qubits, n_elec) + + excitations = sort_excitations(generate_excitations(2, list(range(n_elec)), list(range(n_elec, n_qubits)))) + excitations += sort_excitations(generate_excitations(1, list(range(n_elec)), list(range(n_elec, n_qubits)))) + + excitation_gradient = {tuple(excitation):0.0 for excitation in excitations} + count = 0 + current_energy = mol.e_hf + n_fixed_params = 0 + # Iterating through all excitations; loop breaks if difference in VQE result between excitations is v small + while excitations: + print(f"Iteration {count+1}:") + for excitation in excitations: + _circuit = circuit.copy() + _circuit += givens_excitation_circuit(n_qubits, excitation) + n_parameters = len(_circuit.get_parameters()) + gradient = [round(parameter_shift(_circuit, hamiltonian, index), 5) for index in range(n_parameters)][n_fixed_params:] + # print(f"Energy gradient for {excitation}: {gradient}") + excitation_gradient[tuple(excitation)] = gradient[0] # Gradient magnitude is equal throughout + + # Find the excitation corresponding to the largest gradient, and add it to the circuit + max_grad = max(excitation_gradient, key=lambda x: abs(excitation_gradient.get(x))) + print(f"Excitation with the largest gradient: {max_grad}; Gradient = {excitation_gradient[max_grad]}") + circuit += givens_excitation_circuit(n_qubits, max_grad) + # Remove max_grad from excitations and excitation_data + excitations.pop(excitations.index(list(max_grad))) + del excitation_gradient[max_grad] + + # Run VQE with the updated circuit + vqe = VQE(circuit, hamiltonian) + + circuit_parameters = [param for _tuple in circuit.get_parameters() for param in _tuple] + best, params, extra = vqe.minimize(circuit_parameters, method='BFGS', compile=False) + + n_fixed_params = len(params) + + print(f"\nExact result: {exact_result:.7f}") + print(f" VQE result: {best:.7f}") + energy_diff = best - current_energy + print(f"Difference to previous result: {energy_diff:.7f}") + + if abs(energy_diff) < 1e-3: + print("\nEnergy has converged; exiting while loop") + break + print() + # Update circuit parameters and current energy + circuit.set_parameters(params) + current_energy = best + count += 1 + + + print("\nFinal circuit:") + print(circuit.draw()) + print("\nCircuit statistics:") + print(circuit.summary()) + +Output: + +.. code-block:: output + + Iteration 1: + Excitation with the largest gradient: (2, 3, 6, 7); Gradient = 0.12225 + + Exact result: -7.8770974 + VQE result: -7.8732886 + Difference to previous result: -0.0127499 + + Iteration 2: + Excitation with the largest gradient: (2, 3, 4, 7); Gradient = -0.03485 + + Exact result: -7.8770974 + VQE result: -7.8748417 + Difference to previous result: -0.0015531 + + Iteration 3: + Excitation with the largest gradient: (2, 3, 5, 6); Gradient = 0.03364 + + Exact result: -7.8770974 + VQE result: -7.8762910 + Difference to previous result: -0.0014493 + + Iteration 4: + Excitation with the largest gradient: (0, 1, 4, 5); Gradient = 0.02124 + + Exact result: -7.8770974 + VQE result: -7.8763762 + Difference to previous result: -0.0000853 + + Energy has converged; exiting while loop + + Final circuit: + q0: ─X──────────────────────────────────────────────────────────────────── ... + q1: ─X──────────────────────────────────────────────────────────────────── ... + q2: ─X───o─H─o───RY─o─────RY───X─RY─────o─RY─o─X─H─o─────o─H─o───RY─o───── ... + q3: ─X───|───X───RY─|───X─RY─X─|─RY─X───|─RY─X─|───|─────|───X───RY─|───X─ ... + q4: ─────|──────────|───|────|─|────|───|──────|───|───o─X─────o────|───|─ ... + q5: ─────|──────────|───|────|─|────|───|──────|───|───|───────|────|───|─ ... + q6: ───o─X─────o────|───|────o─o────|───|──────o───X─o─|───────|────|───|─ ... + q7: ───X───H───X────X─H─o───────────o─H─X────────H───X─X───H───X────X─H─o─ ... + + q0: ... ────────────────────────────────────────────────────────────────────── ... + q1: ... ────────────────────────────────────────────────────────────────────── ... + q2: ... RY───X─RY─────o─RY─o─X─H─o─────o─H─o───RY─o─────RY───X─RY─────o─RY─o─X ... + q3: ... RY─X─|─RY─X───|─RY─X─|───|─────|───X───RY─|───X─RY─X─|─RY─X───|─RY─X─| ... + q4: ... ───o─o────|───|──────o───X─o───|──────────|───|────|─|────|───|──────| ... + q5: ... ──────────|───|────────────|─o─X─────o────|───|────o─o────|───|──────o ... + q6: ... ──────────|───|────────────|─X───H───X────X─H─o───────────o─H─X─────── ... + q7: ... ──────────o─H─X────────H───X────────────────────────────────────────── ... + + q0: ... ─────────o─H─o───RY─o─────RY───X─RY─────o─RY─o─X─H─o─── + q1: ... ─────────|───X───RY─|───X─RY─X─|─RY─X───|─RY─X─|───|─── + q2: ... ─H─o─────|──────────|───|────|─|────|───|──────|───|─── + q3: ... ───|─────|──────────|───|────|─|────|───|──────|───|─── + q4: ... ───|───o─X─────o────|───|────o─o────|───|──────o───X─o─ + q5: ... ───X─o─X───H───X────X─H─o───────────o─H─X────────H───X─ + q6: ... ─H───X───────────────────────────────────────────────── + q7: ... ─────────────────────────────────────────────────────── + + Circuit statistics: + Circuit depth = 78 + Total number of gates = 116 + Number of qubits = 8 + Most common gates: + cx: 56 + ry: 32 + h: 24 + x: 4 + +Recall that the full UCCSD circuit for our system had a circuit depth of 1874, with more than 1000 CNOT gates required. +In contrast, the use of a simpler circuit ansatz in conjunction with an adaptive approach allowed us to find a VQE energy that is within chemical accuracy, +while using only 56 CNOT gates and with a final gate depth of only 78. +This is a >20-fold reduction in the gate depth and number of CNOT gates used! .. comment diff --git a/doc/source/tutorials/measurement.rst b/doc/source/tutorials/measurement.rst index 58c18b9..9bc817d 100644 --- a/doc/source/tutorials/measurement.rst +++ b/doc/source/tutorials/measurement.rst @@ -42,9 +42,12 @@ In practice, the expectation value for each of the individual Pauli terms have t This process of obtaining the electronic energy (Hamiltonian expectation value) is still reasonable for a small system. However, the number of Pauli terms in a molecular Hamiltonian scales on the order of :math:`O(N^4)`, where N is the number of qubits. +.. warning:: + + The code block below might take a few minutes to run! + .. code-block:: python - # Warning: This code block might take a few minutes to run from qibochem.driver import Molecule # Build the N2 molecule and get the molecular Hamiltonian From 5a6e45a1ab386a37f1a83ffa164b44e88b6ee902 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Thu, 22 Aug 2024 09:00:55 +0000 Subject: [PATCH 28/32] Remove students' old scripts from examples/ Nicer example almost done in Tutorial documentation --- examples/paper3_caleb.py | 122 ------------------- examples/paper3_conan.py | 109 ----------------- examples/psr_paper2.py | 235 ----------------------------------- examples/psr_paper3.py | 257 --------------------------------------- 4 files changed, 723 deletions(-) delete mode 100644 examples/paper3_caleb.py delete mode 100644 examples/paper3_conan.py delete mode 100644 examples/psr_paper2.py delete mode 100644 examples/psr_paper3.py diff --git a/examples/paper3_caleb.py b/examples/paper3_caleb.py deleted file mode 100644 index 122bc37..0000000 --- a/examples/paper3_caleb.py +++ /dev/null @@ -1,122 +0,0 @@ -""" -A rough Python script attempting to implement the VQE circuit ansatz proposed in Arrazola et al. using Qibo. - -Reference paper: https://doi.org/10.22331/q-2022-06-20-742 - -Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Caleb Seow from -Eunoia Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for -Junior College students. -""" - -import numpy as np -from qibo import Circuit, gates, models - -from qibochem.ansatz.hf_reference import hf_circuit -from qibochem.driver.molecule import Molecule -from qibochem.measurement.expectation import expectation - -mol = Molecule(xyz_file="h2.xyz") -mol.run_pyscf() -n = 4 - -# Single Excitation - -g1 = hf_circuit(mol.nso, 2) - -theta = 0.1 - - -def single_ex(a, b): - g = Circuit(4) - g.add(gates.CNOT(a, b)) - g.add(gates.RY(a, theta)) - g.add(gates.CNOT(b, a)) - g.add(gates.RY(a, -theta)) - g.add(gates.CNOT(b, a)) - g.add(gates.CNOT(a, b)) - return g - - -for i in range(2): - for j in range(2, 4): - if (i + j) % 2 == 0: - g1 += single_ex(i, j) - -print(g1.summary()) -print(g1.draw()) - -theta = 0.1 - - -def double_ex(): - g = Circuit(4) - g.add(gates.CNOT(2, 3)) - g.add(gates.CNOT(0, 2)) - g.add(gates.H(0)) - g.add(gates.H(3)) - g.add(gates.CNOT(0, 1)) # 2 simultaneous CNOTs? - g.add(gates.CNOT(2, 3)) - g.add(gates.RY(0, -theta)) - g.add(gates.RY(1, theta)) - g.add(gates.CNOT(0, 3)) - g.add(gates.H(3)) - g.add(gates.CNOT(3, 1)) - g.add(gates.RY(0, -theta)) - g.add(gates.RY(1, theta)) - g.add(gates.CNOT(2, 1)) - g.add(gates.CNOT(2, 0)) - g.add(gates.RY(0, theta)) - g.add(gates.RY(1, -theta)) - g.add(gates.CNOT(3, 1)) - g.add(gates.H(3)) - g.add(gates.CNOT(0, 3)) - g.add(gates.RY(0, theta)) - g.add(gates.RY(1, -theta)) - g.add(gates.CNOT(0, 1)) - g.add(gates.CNOT(2, 0)) - g.add(gates.H(0)) - g.add(gates.H(3)) - g.add(gates.CNOT(0, 2)) - g.add(gates.CNOT(2, 3)) - return g - - -g1 += double_ex() - -print(g1.draw()) - -hamiltonian = mol.hamiltonian() -# parameter = float(input("theta = ")) -parameter = 0.1 -params = np.zeros(len(g1.get_parameters())) + parameter -g1.set_parameters(params) - -energy = expectation(g1, hamiltonian) -print(f"Expectation: {energy}") - -vqe = models.VQE(g1, hamiltonian) - -best, params, extra = vqe.minimize(params, method="BFGS", compile=False) -# Methods: BFGS, COBYLA, Powell - -print(f"Energy: {best}") -exact_result = hamiltonian.eigenvalues()[0] -print(f"Exact result: {exact_result}") - -hamiltonian = mol.hamiltonian() -# parameter = float(input("theta = ")) -parameter = 0.1 -params = np.zeros(len(g1.get_parameters())) + parameter -g1.set_parameters(params) - -energy = expectation(g1, hamiltonian) -print(f"Expectation: {energy}") - -vqe = models.VQE(g1, hamiltonian) - -best, params, extra = vqe.minimize(params, method="BFGS", compile=False) -# Methods: BFGS, COBYLA, Powell - -print(f"Energy: {best}") - -print(f"Exact result: {exact_result}") diff --git a/examples/paper3_conan.py b/examples/paper3_conan.py deleted file mode 100644 index 4ebcd23..0000000 --- a/examples/paper3_conan.py +++ /dev/null @@ -1,109 +0,0 @@ -""" -A rough Python script attempting to implement the VQE circuit ansatz proposed in Arrazola et al. using Qibo. - -Reference paper: https://doi.org/10.22331/q-2022-06-20-742 - -Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Conan Tan from -National Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for -Junior College students. -""" - -import numpy as np -from qibo import Circuit, gates, models - -from qibochem.ansatz.hf_reference import hf_circuit -from qibochem.driver.molecule import Molecule -from qibochem.measurement.expectation import expectation - -# Excitation function - - -def checkexcitation(elec, orb): - s_excitations = [(i, j) for i in range(elec) for j in range(elec, orb) if (i + j) % 2 == 0] - print(s_excitations) - - d_excitations = [ - (0, 1, k, l) - for k in range(elec, orb) - for l in range(k + 1, orb) - if (1 + k + l) % 2 == 0 and ((k % 2) + (l % 2)) == 1 - ] - print(d_excitations) - return s_excitations, d_excitations - - -# H2 - -# Define molecule and populate -mol = Molecule(xyz_file="h2.xyz") -mol.run_pyscf() -n_qubits = mol.nso -n_electrons = 2 -hamiltonian = mol.hamiltonian() -s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) - - -## Circuit construction -c = hf_circuit(n_qubits, n_electrons) - - -def build(c, singleex, doubleex, x): - for qa, qb in singleex: - sc = Circuit(n_qubits) - sc.add(gates.CNOT(qa, qb)) - sc.add(gates.RY(qa, theta=x / 2)) - sc.add(gates.CNOT(qb, qa)) - sc.add(gates.RY(qa, theta=-x / 2)) - sc.add(gates.CNOT(qb, qa)) - sc.add(gates.CNOT(qa, qb)) - c += sc - # for _i, parameter in enumerate(c.get_parameters()): - # gradient = parameter_shift(c, hamiltonian, parameter_index=_i) - # print(f"Excitation {qa, qb} => Gradient: {gradient}") - # if np.abs(gradient) > 1e-10: - # gradients[(qa, qb)] = gradient - # break - for qa, qb, qc, qd in doubleex: - dc = Circuit(n_qubits) - dc.add(gates.CNOT(qc, qd)) - dc.add(gates.CNOT(qa, qc)) - dc.add(gates.H(qa)) - dc.add(gates.H(qd)) - dc.add(gates.CNOT(qa, qb)) - dc.add(gates.CNOT(qc, qd)) - dc.add(gates.RY(qa, theta=-x / 8)) - dc.add(gates.RY(qb, theta=x / 8)) - dc.add(gates.CNOT(qa, qd)) - dc.add(gates.H(qd)) - dc.add(gates.CNOT(qd, qb)) - dc.add(gates.RY(qa, theta=-x / 8)) - dc.add(gates.RY(qb, theta=x / 8)) - dc.add(gates.CNOT(qc, qb)) - dc.add(gates.CNOT(qc, qa)) - dc.add(gates.RY(qa, theta=x / 8)) - dc.add(gates.RY(qb, theta=-x / 8)) - dc.add(gates.CNOT(qd, qb)) - dc.add(gates.H(qd)) - dc.add(gates.CNOT(qa, qd)) - dc.add(gates.RY(qa, theta=x / 8)) - dc.add(gates.RY(qb, theta=-x / 8)) - dc.add(gates.CNOT(qa, qb)) - dc.add(gates.CNOT(qc, qa)) - dc.add(gates.H(qa)) - dc.add(gates.H(qd)) - dc.add(gates.CNOT(qa, qc)) - dc.add(gates.CNOT(qc, qd)) - c += dc - # for _i, parameter in enumerate(c.get_parameters()): - # gradient = parameter_shift(c, hamiltonian, parameter_index=_i) - # print(f"Excitation {qa, qb, qc, qd} => Gradient: {gradient}") - # if np.abs(gradient) > 1e-10: - # gradients[(qa, qb, qc, qd)] = gradient - # break - return c - - -print(list(enumerate(c.get_parameters()))) -c = build(c, s_excitations, d_excitations, 0.1) -print(c.draw()) -print(c.summary()) diff --git a/examples/psr_paper2.py b/examples/psr_paper2.py deleted file mode 100644 index 37ddbb4..0000000 --- a/examples/psr_paper2.py +++ /dev/null @@ -1,235 +0,0 @@ -""" -This script is a continuation of 'paper2.py'. After implementing the circuit ansatz, this script attempts to examine the -use of the Parameter Shift Rule implemented in Qibo for chemistry applications. In addition, simulations of LiH are used -for some benchmarking. - -Reference paper: https://doi.org/10.1038/s41534-019-0240-1 - -Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Caleb Seow from -Eunoia Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for -Junior College students. -""" - -import math - -import matplotlib.pyplot as plt -import numpy as np -from qibo import Circuit, gates, models -from qibo.derivative import parameter_shift -from qibo.optimizers import optimize -from scipy.optimize import minimize - -from qibochem.ansatz import hf_circuit, ucc_circuit -from qibochem.driver import Molecule -from qibochem.measurement import expectation - -# mol = Molecule(xyz_file="h2.xyz") -h2 = Molecule([("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.4))]) -h2.run_pyscf() - -h2.hf_embedding(active=[1, 2, 5]) - -print(f"# electrons: {h2.n_active_e}") -print(f"# qubits: {h2.n_active_orbs}") - -# circuit = hf_circuit(h2.n_active_orbs, h2.n_active_e) -# print(circuit.draw()) -# print(h2.eps) -# quit() - -hamiltonian = h2.hamiltonian(oei=h2.embed_oei, tei=h2.embed_tei, constant=h2.inactive_energy) - -# Exact_result = hamiltonian.eigenvalues()[0] -# print(f"Exact result: {Exact_result}") -# quit() - -theta, phi = 0.0, 0.0 - - -def special_R(a, theta, phi): - A = Circuit(h2.n_active_orbs) - A.add(gates.RY(a, theta)) - A.add(gates.RZ(a, phi)) - return A - - -def A_gate(a): - A = Circuit(h2.n_active_orbs) - A.add(gates.CNOT(a + 1, a)) - A += special_R(a + 1, theta, phi).invert() - A.add(gates.CNOT(a, a + 1)) - A += special_R(a + 1, theta, phi) - A.add(gates.CNOT(a + 1, a)) - return A - - -# print(circuit.draw()) -# print(A_gate(0).get_parameters()) # 4 gradients - -# initial state -if h2.n_active_e <= math.ceil(h2.n_active_orbs): - Xpositions = [2 * i for i in range(h2.n_active_e)] -else: - Xpositions = [2 * i for i in range(math.ceil(h2.n_active_orbs / 2))] - Xpositions += [2 * j + 1 for j in range(h2.n_active_e - math.ceil(h2.n_active_orbs / 2))] - Xpositions.sort() - -print(Xpositions) -quit() - - -# setting the circuit up -def setup_electrons(g2): - for i in Xpositions: - g2.add(gates.X(i)) - - -Xpositions = [0, 2] - -print("Xpositions: ", Xpositions) - -# Xpositions is like a starting point -no_of_gates = 0 -var_list = [] -master_list = [] -ref_list = [] -popping_list = [] - -# initialising -for i in range(h2.n_active_orbs): - if i + 1 not in Xpositions: - if i != h2.n_active_orbs - 1: - ref_list.append((i, i + 1)) - master_list.append((i, i + 1)) - no_of_gates += 1 - -print(f"initial master list: {master_list}") - - -# Adding 1 layer based on the prev layer -while no_of_gates < math.comb(h2.n_active_orbs, h2.n_active_e): - - for i, j in ref_list: - var_list.append((i - 1, i)) - var_list.append((j, j + 1)) - - var_list.append((0, 0)) # dummy tuple that won't appear naturally - - for i in range(len(var_list) - 1): # 3 conditions to pop - if var_list[i][0] < 0: - popping_list.append(i) - elif var_list[i][1] > h2.n_active_orbs - 1: - popping_list.append(i) - elif var_list[i] == var_list[i + 1]: - popping_list.append(i + 1) - - popping_list.reverse() # so that index won't change after every pop - - for i in popping_list: - var_list.pop(i) - - var_list.pop() # removing the (0,0) tuple - - for i in var_list: - master_list.append(i) - no_of_gates += 1 - - ref_list.clear() - for i in var_list: - ref_list.append(i) - - popping_list.clear() - var_list.clear() - -print(f"final master list: {master_list}") - -A_gate_indexes = master_list[: math.comb(h2.n_active_orbs, h2.n_active_e)] # shouldnt be needed but jic?? - -print(f"A gate indexes: {A_gate_indexes}") - -g2 = Circuit(h2.n_active_orbs) - -setup_electrons(g2) - -for i, j in A_gate_indexes: - g2 += A_gate(i) - -print(g2.draw()) - -t, p = 0.1, 0.1 # theta, phi -params = [] - -# qn: for the dagger do we want the changed params too? -# A_gate_param = [-p-np.pi,-t-np.pi/2,t+np.pi/2,p+np.pi] -# A_gate_param = [p+np.pi,t+np.pi/2,t+np.pi/2,p+np.pi] -A_gate_param = [-t, -p, p, t] - -for i in range(len(A_gate_indexes)): - params += A_gate_param - -g2.set_parameters(params) - -# gradlist = [] -# -# gradients = {} -# -# print(len(g2.get_parameters())) -# -# thetalist = [f"theta {i}" for i in range(len(g2.get_parameters()))] -# -# for _i, parameter in enumerate(g2.get_parameters()): -# gradient = parameter_shift(g2, hamiltonian, parameter_index=_i) -# gradlist.append(gradient) -# -# for i,j in enumerate(thetalist): -# if np.abs(gradlist[i]) > 1e-10: -# gradients[j] = gradient - -# print(gradients) - -energy = expectation(g2, hamiltonian) -print(f"Expectation: {energy}") - -# vqe = models.VQE(g2, hamiltonian) - -# Methods: BFGS, COBYLA, Powell -# best_BFGS, params_BFGS, extra_BFGS = vqe.minimize(params, method='BFGS', compile=False) -# print(f"Energy BFGS: {best_BFGS}") -# -# best_COBYLA, params_COBYLA, extra_COBYLA = vqe.minimize(params, method='COBYLA', compile=False) -# print(f"Energy COBYLA: {best_COBYLA}") -# -# best_Powell, params_Powell, extra_Powell = vqe.minimize(params, method='Powell', compile=False) -# print(f"Energy Powell: {best_Powell}") - -# Exact_result = hamiltonian.eigenvalues()[0] -# print(f"Exact result: {Exact_result}") - - -def no_restriction_electronic_energy(parameters): - g2.set_parameters(parameters) - return expectation(g2, hamiltonian) - - -# Combine the gradients of each parameterised gate into a single array (vector): -def gradient_no_restriction(_thetas): - g2.set_parameters(_thetas) - return np.array([parameter_shift(g2, hamiltonian, parameter_index=_i) for _i in range(len(_thetas))]) - - -# Start from zero -thetas = np.zeros(len(g2.get_parameters())) - -best, params, extra = optimize(no_restriction_electronic_energy, thetas, method="SLSQP", jac=gradient_no_restriction) -# best, params, extra = vqe.minimize(thetas, method='BFGS', jac=ucc_gradient_no_restriction) - -print(f"VQE result: {best}") -print(f"Optimized parameter: {params}") -print(f"Number of steps: {extra.nfev}") -print() - -# COBYLA -best, params, extra = optimize(no_restriction_electronic_energy, thetas, method="SLSQP") -print(f"VQE result: {best}") -print(f"Optimized parameter: {params}") -print(f"Number of steps: {extra.nfev}") diff --git a/examples/psr_paper3.py b/examples/psr_paper3.py deleted file mode 100644 index 4ae8f8c..0000000 --- a/examples/psr_paper3.py +++ /dev/null @@ -1,257 +0,0 @@ -""" -This script is a continuation of 'paper3_conan.py'. After implementing the circuit ansatz, this script attempts to -examine the use of the Parameter Shift Rule implemented in Qibo for chemistry applications. In addition, simulations of -LiH are used for some benchmarking. - -Reference paper: https://doi.org/10.22331/q-2022-06-20-742 - -Acknowledgements: The original draft of this code, in the form of a Jupyter notebook, was prepared by Conan Tan from -National Junior College, who was attached to IHPC in December 2023 under the A*STAR Research Attachment Programme for -Junior College students. -""" - -import numpy as np -from qibo import Circuit, gates -from qibo.derivative import parameter_shift -from qibo.optimizers import optimize - -from qibochem.ansatz.hf_reference import hf_circuit -from qibochem.driver.molecule import Molecule -from qibochem.measurement.expectation import expectation - - -# Excitation function -def checkexcitation(elec, orb): - s_excitations = [(i, j) for i in range(elec) for j in range(elec, orb) if (i + j) % 2 == 0] - print(s_excitations) - - d_excitations = [ - (0, 1, k, l) - for k in range(elec, orb) - for l in range(k + 1, orb) - if (1 + k + l) % 2 == 0 and ((k % 2) + (l % 2)) == 1 - ] - print(d_excitations) - return s_excitations, d_excitations - - -# H2 -# Define molecule and populate -mol = Molecule(xyz_file="h2.xyz") -mol.run_pyscf() -n_qubits = mol.nso -n_electrons = mol.nelec -hamiltonian = mol.hamiltonian() -s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) -gradients = {} - -## Circuit construction -c = hf_circuit(n_qubits, n_electrons) - - -def build(c, singleex, doubleex, x): - for qa, qb in singleex: - sc = Circuit(n_qubits) - sc.add(gates.CNOT(qa, qb)) - sc.add(gates.RY(qa, theta=x / 2)) - sc.add(gates.CNOT(qb, qa)) - sc.add(gates.RY(qa, theta=-x / 2)) - sc.add(gates.CNOT(qb, qa)) - sc.add(gates.CNOT(qa, qb)) - c += sc - for _i, parameter in enumerate(c.get_parameters()): - gradient = parameter_shift(c, hamiltonian, parameter_index=_i) - print(f"Excitation {qa, qb} => Gradient: {gradient}") - if np.abs(gradient) > 1e-10: - gradients[(qa, qb)] = gradient - # break - - # for qa, qb, qc, qd in doubleex: - # dc = Circuit(n_qubits) - # dc.add(gates.CNOT(qc, qd)) - # dc.add(gates.CNOT(qa, qc)) - # dc.add(gates.H(qa)) - # dc.add(gates.H(qd)) - # dc.add(gates.CNOT(qa, qb)) - # dc.add(gates.CNOT(qc, qd)) - # dc.add(gates.RY(qa, theta=-x/8)) - # dc.add(gates.RY(qb, theta=x/8)) - # dc.add(gates.CNOT(qa, qd)) - # dc.add(gates.H(qd)) - # dc.add(gates.CNOT(qd, qb)) - # dc.add(gates.RY(qa, theta=-x/8)) - # dc.add(gates.RY(qb, theta=x/8)) - # dc.add(gates.CNOT(qc, qb)) - # dc.add(gates.CNOT(qc, qa)) - # dc.add(gates.RY(qa, theta=x/8)) - # dc.add(gates.RY(qb, theta=-x/8)) - # dc.add(gates.CNOT(qd, qb)) - # dc.add(gates.H(qd)) - # dc.add(gates.CNOT(qa, qd)) - # dc.add(gates.RY(qa, theta=x/8)) - # dc.add(gates.RY(qb, theta=-x/8)) - # dc.add(gates.CNOT(qa, qb)) - # dc.add(gates.CNOT(qc, qa)) - # dc.add(gates.H(qa)) - # dc.add(gates.H(qd)) - # dc.add(gates.CNOT(qa, qc)) - # dc.add(gates.CNOT(qc, qd)) - # c += dc - # for _i, parameter in enumerate(c.get_parameters()): - # gradient = parameter_shift(c, hamiltonian, parameter_index=_i) - # print(f"Excitation {qa, qb, qc, qd} => Gradient: {gradient}") - # if np.abs(gradient) > 1e-10: - # gradients[(qa, qb, qc, qd)] = gradient - # break - return c - - -# checkexcitation(sum(n_electrons), n_qubits) -print(list(enumerate(c.get_parameters()))) -c = build(c, s_excitations, d_excitations, 0.0) -print(c.draw()) -print(c.summary()) - -for excitation in sorted(gradients, key=lambda x: np.abs(gradients[x]), reverse=True): - print(f"Excitation {excitation} => Gradient: {gradients[excitation]}") - -params = [] -s_coeffs = [1 / 2, -1 / 2] -d_coeffs = [-1 / 8, 1 / 8, -1 / 8, 1 / 8, 1 / 8, -1 / 8, 1 / 8, -1 / 8] - -for i in s_excitations: - params += s_coeffs - -for i in d_excitations: - params += d_coeffs - -## Restricted -x = [] - - -def electronic_energy(x): - x *= params - c.set_parameters(x) - return expectation(c, hamiltonian) - - -y = [] - - -def gradient_restriction(y): - y *= params - c.set_parameters(y) - return np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))]) - - -print(np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))])) - - -## Non-Restricted -def nr_electronic_energy(parameter): - c.set_parameters(parameter) - return expectation(c, hamiltonian) - - -y = [] - - -def gradient_no_restriction(y): - return np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))]) - - -print(np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))])) - - -## Energy Calculation -methods = ["SLSQP", "BFGS", "COBYLA", "Powell"] - -for m in methods: - theta = np.zeros(len(c.get_parameters())) - if m in ["BFGS", "SLSQP"]: - # Gradients for BFGS and SLSQP - best, params, extra = optimize(electronic_energy, theta, method=m, jac=gradient_restriction) - else: - best, params, extra = optimize(electronic_energy, theta, method=m) - - print(f"Energy {m}: {best}") - print(f"Optimized parameter: {params}") - print(f"Number of steps: {extra.nfev}") - -print(f"Exact result: {hamiltonian.eigenvalues()[0]}") - - -# LiH -mol = Molecule(xyz_file="lih.xyz") -mol.run_pyscf() -mol.hf_embedding(active=[1, 2, 5]) -hamiltonian = mol.hamiltonian(oei=mol.embed_oei, tei=mol.embed_tei, constant=mol.inactive_energy) -n_qubits = mol.n_active_orbs -n_electrons = mol.n_active_e -s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) -gradients = {} - -## Circuit construction - -c = hf_circuit(n_qubits, n_electrons) -s_excitations, d_excitations = checkexcitation(n_electrons, n_qubits) -c = build(c, s_excitations, d_excitations, 0.1) -print(c.draw()) -print(c.summary()) - -print(f"\nInitial number of excitation: {len(s_excitations + d_excitations)}") -print(f"Final number of excitations: {len(gradients)}") - -for excitation in sorted(gradients, key=lambda x: np.abs(gradients[x]), reverse=True): - print(f"Excitation {excitation} => Gradient: {gradients[excitation]}") - -## Calculation of Hamiltonian - -param_range = np.arange(*map(float, input("Enter space-seperated parameter range, and step increment = ").split())) -for x in param_range: - print(f"Theta: {round(x, 5)} => Energy expectation value: {electronic_energy(x)}") - -## Restricted -x = [] - - -def electronic_energy(x): - x *= params - c.set_parameters(x) - return expectation(c, hamiltonian) - - -# param_range = np.arange(*map(float, input("Enter space-seperated parameter range, and step increment = ").split())) -# for x in param_range: -# print(f"Theta: {round(x, 5)} => Energy expectation value: {electronic_energy(x)}") - -# print(params) - -y = [] - - -def gradient_restriction(y): - y *= params - c.set_parameters(y) - return np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))]) - - -print(np.array([parameter_shift(c, hamiltonian, parameter_index=_i) for _i in range(len(y))])) - -## Energy Calculation - -methods = ["SLSQP", "BFGS", "COBYLA", "Powell"] - -for m in methods: - theta = np.zeros(len(c.get_parameters())) - if m in ["BFGS", "SLSQP"]: - # Gradients for BFGS and SLSQP - best, params, extra = optimize(electronic_energy, theta, method=m, jac=gradient_restriction) - else: - best, params, extra = optimize(electronic_energy, theta, method=m) - - print(f"Energy {m}: {best}") - print(f"Optimized parameter: {params}") - print(f"Number of steps: {extra.nfev}") - -print(f"Exact result: {hamiltonian.eigenvalues()[0]}") From 4e81792655c6d09fcbf879d353b86f915159d117 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Fri, 13 Sep 2024 02:52:17 +0000 Subject: [PATCH 29/32] Add references and clean up tutorial --- doc/source/tutorials/adaptive.rst | 241 +----------------------------- 1 file changed, 5 insertions(+), 236 deletions(-) diff --git a/doc/source/tutorials/adaptive.rst b/doc/source/tutorials/adaptive.rst index 9dbf19e..f3189f2 100644 --- a/doc/source/tutorials/adaptive.rst +++ b/doc/source/tutorials/adaptive.rst @@ -32,12 +32,12 @@ Output: As shown in the above code block, the full UCCSD circuit for a simplified LiH/STO-3G system has a circuit depth of 1874 (!), with more than 1000 CNOT gates required! Hence, there is still a need to further reduce and simplify the circuit ansatzes used for running a VQE simulation. -Other than designing shorter and more efficient circuit ansatzes, one alternative approach is through the use of energy gradients - for instance, through the Parameter Shift Rule on hardware - to filter and reduce the number of fermionic excitations in a circuit ansatz. (REFS!!!) +Other than designing shorter and more efficient circuit ansatzes, one alternative approach is through the use of energy gradients - for instance, through the Parameter Shift Rule on hardware - to filter and reduce the number of fermionic excitations in a circuit ansatz. [#f1]_ [#f2]_ This is known as an adaptive method, in the sense that the quantum gates used to construct the circuit ansatz, as well as its actual structure and arrangement is not fixed, and varies depending on the molecular system under study. For example, in a H2/STO-3G system mapped with the Jordan-Wigner transformation, there are three possible spin-allowed fermionic excitations: two single excitations (``[0, 2]``, ``[1, 3]``) and one double excitation (``[0, 1, 2, 3]``). -The full UCCSD circuit for this system has been shown in an earlier example (ref), and it requires 64 CNOT gates for this simple molecular system. +The full UCCSD circuit for this system has been shown in an earlier :ref:`example `, and it requires 64 CNOT gates for this simple molecular system. Let's look at the gradients of each of the individual fermionic excitations: @@ -108,7 +108,7 @@ Output: We managed to find the exact result by applying only the doubles excitation! Next, let's look at the potential savings for the simplified LiH/STO-3G system. -To reduce the circuit depth further, we will use the more modern ansatz, the Givens excitation circuit from Arrazola et al., instead of the UCC ansatz. +To reduce the circuit depth further, we will use the more modern ansatz, the Givens excitation circuit from Arrazola et al., [#f1]_ instead of the UCC ansatz. As was done in the above example, we will start with a HF circuit, then find the gradients for each circuit ansatz corresponding to a fermionic excitation. After that, the excitation with the largest absolute value of the gradient will be added to the initial circuit, followed by a VQE simulation. @@ -359,239 +359,8 @@ while using only 56 CNOT gates and with a final gate depth of only 78. This is a >20-fold reduction in the gate depth and number of CNOT gates used! -.. comment - OLD STUFF - --------- - - A quantum circuit comprising parameterized gates (`e.g.` :math:`RX(\theta)`, :math:`RY(\theta)` and :math:`RZ(\theta)`), - represents a unitary transformation :math:`U(\theta)` that transforms some initial quantum state into a parametrized ansatz state :math:`|\psi(\theta)\rangle`. - - Examples of some ansatzes available in Qibochem are described in the subsections below. - - Hardware Efficient Ansatz - ------------------------- - - Qibochem provides a hardware efficient ansatz that simply consists of a layer of single-qubit rotation gates followed by a layer of two-qubit gates that entangle the qubits. - For the H\ :sub:`2` case discussed in previous sections, a possible hardware efficient circuit ansatz can be constructed as such: - - .. image:: qibochem_doc_ansatz_hardware-efficient.svg - - .. code-block:: python - - from qibochem.ansatz import he_circuit - - nqubits = 4 - nlayers = 1 - - circuit = he_circuit(nqubits, nlayers) - print(circuit.draw()) - - .. code-block:: output - - q0: ─RY─RZ─o─────Z─ - q1: ─RY─RZ─Z─o───|─ - q2: ─RY─RZ───Z─o─|─ - q3: ─RY─RZ─────Z─o─ - - The energy of the state generated from the hardware efficient ansatz for the fermionic two-body Hamiltonian can then be estimated, using state vectors or samples. - - The following example demonstrates how the energy of the H2 molecule is affected with respect to the rotational parameters: - - .. code-block:: python - - import numpy as np - - from qibochem.driver import Molecule - from qibochem.measurement.expectation import expectation - 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() - hamiltonian = mol.hamiltonian() - - # Define and build the HEA - nlayers = 1 - nqubits = mol.nso - ntheta = 2 * nqubits * nlayers - circuit = he_circuit(nqubits, nlayers) - - print("Energy expectation values for thetas: ") - print("-----------------------------") - print("| theta | Electronic energy |") - print("|---------------------------|") - thetas = [-0.2, 0.0, 0.2] - for theta in thetas: - params = np.full(ntheta, theta) - circuit.set_parameters(params) - electronic_energy = expectation(circuit, hamiltonian) - print(f"| {theta:5.1f} | {electronic_energy:^18.12f}|") - print("-----------------------------") - - - .. code-block:: output - - converged SCF energy = -1.11628373627429 - - Energy expectation values for thetas: - ----------------------------- - | theta | Electronic energy | - |---------------------------| - | -0.2 | 0.673325849299 | - | 0.0 | 0.707418334474 | - | 0.2 | 0.673325849299 | - ----------------------------- - - - .. _UCC Ansatz: - - Unitary Coupled Cluster Ansatz - ------------------------------ - - The Unitary Coupled Cluster (UCC) ansatz [#f1]_ [#f2]_ [#f3]_ is a variant of the popular gold standard Coupled Cluster ansatz [#f4]_ of quantum chemistry. - The UCC wave function is a parameterized unitary transformation of a reference wave function :math:`\psi_{\mathrm{ref}}`, of which a common choice is the Hartree-Fock wave function. - - .. math:: - - \begin{align*} - |\psi_{\mathrm{UCC}}\rangle &= U(\theta)|\psi_{\mathrm{ref}}\rangle \\ - &= e^{\hat{T}(\theta) - \hat{T}^\dagger(\theta)}|\psi_{\mathrm{ref}}\rangle - \end{align*} - - - Similar to the process for the molecular Hamiltonian, the fermionic excitation operators :math:`\hat{T}` and :math:`\hat{T}^\dagger` are mapped using e.g. Jordan-Wigner mapping into Pauli operators. - This is typically followed by a Suzuki-Trotter decomposition of the exponentials of these Pauli operators, which allows the UCC ansatz to be implemented on quantum computers. [#f5]_ - - An example of how to build a UCC doubles circuit ansatz for the :math:`H_2` molecule is given as: - - .. code-block:: python - - from qibochem.driver.molecule import Molecule - 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() - hamiltonian = mol.hamiltonian() - - # Set parameters for the rest of the experiment - n_qubits = mol.nso - n_electrons = mol.nelec - - # Build UCCD circuit - circuit = hf_circuit(n_qubits, n_electrons) # Start with HF circuit - circuit += ucc_circuit(n_qubits, [0, 1, 2, 3]) # Then add the double excitation circuit ansatz - - print(circuit.draw()) - - .. code-block:: output - - q0: ─X──H─────X─RZ─X─────H──RX─────X─RZ─X─────RX─RX─────X─RZ─X─────RX─H─── ... - q1: ─X──H───X─o────o─X───H──RX───X─o────o─X───RX─H────X─o────o─X───H──RX── ... - q2: ─RX───X─o────────o─X─RX─RX─X─o────────o─X─RX─H──X─o────────o─X─H──H──X ... - q3: ─H────o────────────o─H──H──o────────────o─H──H──o────────────o─H──H──o ... - - q0: ... ───X─RZ─X─────H──RX─────X─RZ─X─────RX─H──────X─RZ─X─────H──H──────X─RZ ... - q1: ... ─X─o────o─X───RX─H────X─o────o─X───H──RX───X─o────o─X───RX─H────X─o─── ... - q2: ... ─o────────o─X─H──RX─X─o────────o─X─RX─RX─X─o────────o─X─RX─H──X─o───── ... - q3: ... ────────────o─H──RX─o────────────o─RX─RX─o────────────o─RX─RX─o─────── ... - - q0: ... ─X─────H──RX─────X─RZ─X─────RX─ - q1: ... ─o─X───H──RX───X─o────o─X───RX─ - q2: ... ───o─X─H──H──X─o────────o─X─H── - q3: ... ─────o─RX─RX─o────────────o─RX─ - - - .. - _Basis rotation ansatz - - Basis rotation ansatz - --------------------- - - The starting points for contemporary quantum chemistry methods are often those based on the mean field approximation within a (finite) molecular orbital basis, i.e. the Hartree-Fock method. The electronic energy is calculated as the mean value of the electronic Hamiltonian :math:`\hat{H}_{\mathrm{elec}}` acting on a normalized single Slater determinant function :math:`\psi` [#f6]_ - - .. math:: - - \begin{align*} - E[\psi] &= \langle \psi | \hat{H}_{\mathrm{elec}} |\psi \rangle \\ - &= \sum_i^{N_f} \langle \phi_i |\hat{h}|\phi_i \rangle + \frac{1}{2} \sum_{i,j}^{N_f} - \langle \phi_i\phi_j||\phi_i\phi_j \rangle - \end{align*} - - The orthonormal molecular orbitals :math:`\phi` are optimized by a direct minimization of the energy functional with respect to parameters :math:`\kappa` that parameterize the unitary rotations of the orbital basis. Qibochem's implementation uses the QR decomposition of the unitary matrix as employed by Clements et al., [#f7]_ which results in a rectangular gate layout of `Givens rotation gates `_ that yield linear CNOT gate depth when decomposed. - - - .. code-block:: python - - import numpy as np - from qibochem.driver.molecule import Molecule - from qibochem.ansatz import basis_rotation, ucc - from qibo import Circuit, gates, models - - def basis_rotation_circuit(mol, parameters=0.0): - - nqubits = mol.nso - occ = range(0, mol.nelec) - vir = range(mol.nelec, mol.nso) - - U, kappa = basis_rotation.unitary(occ, vir, parameters=parameters) - gate_angles, final_U = basis_rotation.givens_qr_decompose(U) - gate_layout = basis_rotation.basis_rotation_layout(nqubits) - gate_list, ordered_angles = basis_rotation.basis_rotation_gates(gate_layout, gate_angles, kappa) - - circuit = Circuit(nqubits) - for _i in range(mol.nelec): - circuit.add(gates.X(_i)) - circuit.add(gate_list) - - return circuit, gate_angles - - h3p = Molecule([('H', (0.0000, 0.0000, 0.0000)), - ('H', (0.0000, 0.0000, 0.8000)), - ('H', (0.0000, 0.0000, 1.6000))], - charge=1, multiplicity=1) - h3p.run_pyscf(max_scf_cycles=1) - - e_init = h3p.e_hf - h3p_sym_ham = h3p.hamiltonian("sym", h3p.oei, h3p.tei, 0.0, "jw") - - hf_circuit, qubit_parameters = basis_rotation_circuit(h3p, parameters=0.1) - - print(hf_circuit.draw()) - - vqe = models.VQE(hf_circuit, h3p_sym_ham) - res = vqe.minimize(qubit_parameters) - - print('energy of initial guess: ', e_init) - print('energy after vqe : ', res[0]) - - .. code-block:: output - - q0: ─X─G─────────G─────────G───────── - q1: ─X─G─────G───G─────G───G─────G─── - q2: ─────G───G─────G───G─────G───G─── - q3: ─────G─────G───G─────G───G─────G─ - q4: ───────G───G─────G───G─────G───G─ - q5: ───────G─────────G─────────G───── - basis rotation: using uniform value of 0.1 for each parameter value - energy of initial guess: -1.1977713400022736 - energy after vqe : -1.2024564133305427 - - - - - - .. rubric:: References - .. [#f1] Kutzelnigg, W. (1977). 'Pair Correlation Theories', in Schaefer, H.F. (eds) Methods of Electronic Structure Theory. Modern Theoretical Chemistry, vol 3. Springer, Boston, MA. - - .. [#f2] Whitfield, J. D. et al., 'Simulation of Electronic Structure Hamiltonians using Quantum Computers', Mol. Phys. 109 (2011) 735. - - .. [#f3] Anand. A. et al., 'A Quantum Computing view on Unitary Coupled Cluster Theory', Chem. Soc. Rev. 51 (2022) 1659. - - .. [#f4] Crawford, T. D. et al., 'An Introduction to Coupled Cluster Theory for Computational Chemists', in Reviews in Computational Chemistry 14 (2007) 33. - - .. [#f5] Barkoutsos, P. K. et al., 'Quantum algorithms for electronic structure calculations: Particle-hole Hamiltonian and optimized wave-function expansions', Phys. Rev. A 98 (2018) 022322. - - .. [#f6] Piela, L. (2007). 'Ideas of Quantum Chemistry'. Elsevier B. V., the Netherlands. + .. [#f1] Arrazola, J. M. et al., 'Universal quantum circuits for quantum chemistry', Quantum, 6, (2022), 742. - .. [#f7] Clements, W. R. et al., 'Optimal Design for Universal Multiport Interferometers', Optica 3 (2016) 1460. + .. [#f2] Schuld M. et al., 'Evaluating analytic gradients on quantum hardware', Phys. Rev. A, 99, (2019), 032331. From 6301dd963f6b323c0b8f7ccdc3d8146d8992c607 Mon Sep 17 00:00:00 2001 From: Wong Zi Cheng <70616433+chmwzc@users.noreply.github.com> Date: Fri, 13 Sep 2024 03:00:50 +0000 Subject: [PATCH 30/32] Remove unused file --- examples/h2.xyz | 4 ---- 1 file changed, 4 deletions(-) delete mode 100644 examples/h2.xyz diff --git a/examples/h2.xyz b/examples/h2.xyz deleted file mode 100644 index 6a4e902..0000000 --- a/examples/h2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 - -H 0.0 0.0 0.0 -H 0.0 0.0 0.7 From d46d5c6e88f26f60db39d083e0b5a5907a2b8d99 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 23 Dec 2024 23:38:17 +0000 Subject: [PATCH 31/32] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/asottile/pyupgrade: v3.19.0 → v3.19.1](https://github.com/asottile/pyupgrade/compare/v3.19.0...v3.19.1) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index e6b6ad9..d4262ba 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -20,6 +20,6 @@ repos: - id: isort args: ["--profile", "black"] - repo: https://github.com/asottile/pyupgrade - rev: v3.19.0 + rev: v3.19.1 hooks: - id: pyupgrade From 9fc04fe79dd06a0455561f3dc29df77f260d22b6 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 24 Dec 2024 03:20:41 +0000 Subject: [PATCH 32/32] Bump jinja2 from 3.1.4 to 3.1.5 Bumps [jinja2](https://github.com/pallets/jinja) from 3.1.4 to 3.1.5. - [Release notes](https://github.com/pallets/jinja/releases) - [Changelog](https://github.com/pallets/jinja/blob/main/CHANGES.rst) - [Commits](https://github.com/pallets/jinja/compare/3.1.4...3.1.5) --- updated-dependencies: - dependency-name: jinja2 dependency-type: indirect ... Signed-off-by: dependabot[bot] --- poetry.lock | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/poetry.lock b/poetry.lock index 76e69b1..66c05bb 100644 --- a/poetry.lock +++ b/poetry.lock @@ -1,4 +1,4 @@ -# This file is automatically @generated by Poetry 1.8.3 and should not be changed by hand. +# This file is automatically @generated by Poetry 1.8.5 and should not be changed by hand. [[package]] name = "alabaster" @@ -937,13 +937,13 @@ testing = ["Django", "attrs", "colorama", "docopt", "pytest (<7.0.0)"] [[package]] name = "jinja2" -version = "3.1.4" +version = "3.1.5" description = "A very fast and expressive template engine." optional = false python-versions = ">=3.7" files = [ - {file = "jinja2-3.1.4-py3-none-any.whl", hash = "sha256:bc5dd2abb727a5319567b7a813e6a2e7318c39f4f487cfe6c89c6f9c7d25197d"}, - {file = "jinja2-3.1.4.tar.gz", hash = "sha256:4a3aee7acbbe7303aede8e9648d13b8bf88a429282aa6122a993f0ac800cb369"}, + {file = "jinja2-3.1.5-py3-none-any.whl", hash = "sha256:aba0f4dc9ed8013c424088f68a5c226f7d6097ed89b246d7749c2ec4175c6adb"}, + {file = "jinja2-3.1.5.tar.gz", hash = "sha256:8fefff8dc3034e27bb80d67c671eb8a9bc424c0ef4c0826edbff304cceff43bb"}, ] [package.dependencies]