diff --git a/.github/workflows/github_actions_automated_testing.yml b/.github/workflows/github_actions_automated_testing.yml index ceade1f28..0635cd5fa 100755 --- a/.github/workflows/github_actions_automated_testing.yml +++ b/.github/workflows/github_actions_automated_testing.yml @@ -40,7 +40,7 @@ jobs: pip install qiskit==0.33.1 # Due to strange behaviour of noise model pip install qulacs pip install amazon-braket-sdk - pip install cirq==0.12.0 + pip install cirq pip install projectq if: always() diff --git a/tangelo/algorithms/variational/tests/test_vqe_solver.py b/tangelo/algorithms/variational/tests/test_vqe_solver.py index c3aaedc09..e910be737 100644 --- a/tangelo/algorithms/variational/tests/test_vqe_solver.py +++ b/tangelo/algorithms/variational/tests/test_vqe_solver.py @@ -12,6 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. +import re import unittest import numpy as np @@ -137,6 +138,31 @@ def test_simulate_qcc_h2(self): energy = vqe_solver.simulate() self.assertAlmostEqual(energy, -1.137270, delta=1e-4) + def test_simulate_vsqs_h2(self): + """Run VQE on H2 molecule, with vsqs ansatz, JW qubit mapping, exact simulator for both molecule input and + qubit_hamiltonian/hini/reference_state input + """ + vqe_options = {"molecule": mol_H2_sto3g, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw", + "verbose": True, "ansatz_options": {"intervals": 3, "time": 3}} + vqe_solver = VQESolver(vqe_options) + vqe_solver.build() + + energy = vqe_solver.simulate() + self.assertAlmostEqual(energy, -1.137270, delta=1e-4) + + qubit_hamiltonian = vqe_solver.qubit_hamiltonian + hini = vqe_solver.ansatz.hini + reference_state = vqe_solver.ansatz.prepare_reference_state() + + vqe_options = {"molecule": None, "qubit_hamiltonian": qubit_hamiltonian, "ansatz": BuiltInAnsatze.VSQS, "qubit_mapping": "jw", + "ansatz_options": {"intervals": 3, "time": 3, "qubit_hamiltonian": qubit_hamiltonian, + "hini": hini, "reference_state": reference_state}} + vqe_solver = VQESolver(vqe_options) + vqe_solver.build() + + energy = vqe_solver.simulate() + self.assertAlmostEqual(energy, -1.137270, delta=1e-4) + def test_simulate_h2_qiskit(self): """Run VQE on H2 molecule, with UCCSD ansatz, JW qubit mapping, initial parameters, exact qiskit simulator. diff --git a/tangelo/algorithms/variational/vqe_solver.py b/tangelo/algorithms/variational/vqe_solver.py index 639e3bf3c..e5fe56729 100644 --- a/tangelo/algorithms/variational/vqe_solver.py +++ b/tangelo/algorithms/variational/vqe_solver.py @@ -35,6 +35,7 @@ from tangelo.toolboxes.ansatz_generator.upccgsd import UpCCGSD from tangelo.toolboxes.ansatz_generator.qmf import QMF from tangelo.toolboxes.ansatz_generator.qcc import QCC +from tangelo.toolboxes.ansatz_generator.vsqs import VSQS from tangelo.toolboxes.ansatz_generator.variational_circuit import VariationalCircuitAnsatz from tangelo.toolboxes.ansatz_generator.penalty_terms import combined_penalty from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies @@ -50,6 +51,7 @@ class BuiltInAnsatze(Enum): UpCCGSD = 4 QMF = 5 QCC = 6 + VSQS = 7 class VQESolver: @@ -186,13 +188,20 @@ def build(self): self.ansatz = QMF(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options) elif self.ansatz == BuiltInAnsatze.QCC: self.ansatz = QCC(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options) + elif self.ansatz == BuiltInAnsatze.VSQS: + self.ansatz = VSQS(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options) else: raise ValueError(f"Unsupported ansatz. Built-in ansatze:\n\t{self.builtin_ansatze}") elif not isinstance(self.ansatz, Ansatz): raise TypeError(f"Invalid ansatz dataype. Expecting instance of Ansatz class, or one of built-in options:\n\t{self.builtin_ansatze}") # Building with a qubit Hamiltonian. - elif not isinstance(self.ansatz, Ansatz): + elif not isinstance(self.ansatz, Ansatz) or self.ansatz not in [BuiltInAnsatze.HEA, BuiltInAnsatze.VSQS]: + if self.ansatz == BuiltInAnsatze.HEA: + self.ansatz = HEA(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options) + elif self.ansatz == BuiltInAnsatze.VSQS: + self.ansatz = VSQS(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options) + else: raise TypeError(f"Invalid ansatz dataype. Expecting a custom Ansatz (Ansatz class).") # Set ansatz initial parameters (default or use input), build corresponding ansatz circuit diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py index a9a2a4c16..68cb31630 100644 --- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py +++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py @@ -83,7 +83,7 @@ def exp_pauliword_to_gates(pauli_word, coef, variational=True, control=None): def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=False, trotter_order=1, control=None, - return_phase=False): + return_phase=False, pauli_order=None): """Generate the exponentiation of a qubit operator in first- or second-order Trotterized form. The algorithm is described in Whitfield 2010 https://arxiv.org/pdf/1001.3855.pdf @@ -94,12 +94,19 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals variational (bool) : Whether the coefficients are variational trotter_order (int): order of trotter approximation, only 1 or 2 are supported. return_phase (bool): Return the global-phase generated + pauli_order (list): The desired pauli_word order defined as a list with elements (pauli_word, coeff) + with corresponding dictionary elements pauli_word: coeff in QubitOperator terms.items() Returns: Circuit: circuit corresponding to exponentiation of qubit operator phase : The global phase of the time evolution if return_phase=True else not included """ - pauli_words = list(qubit_op.terms.items()) + if pauli_order is None: + pauli_words = list(qubit_op.terms.items()) + elif isinstance(pauli_order, list): + pauli_words = deepcopy(pauli_order) + else: + raise ValueError("ordered terms must be a list with elements (keys, values) of qubit_op.terms.items()") if trotter_order > 2: raise ValueError(f"Trotter order of >2 is not supported currently in Tangelo.") @@ -118,7 +125,7 @@ def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=Fals phase = 1. exp_pauli_word_gates = list() for i in range(trotter_order): - if i == 0: + if i == 1: pauli_words.reverse() for pauli_word, coef in pauli_words: if pauli_word: # identity terms do not contribute to evolution outside of a phase diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py new file mode 100644 index 000000000..d53e6bf02 --- /dev/null +++ b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py @@ -0,0 +1,116 @@ +# Copyright 2021 Good Chemistry Company. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import unittest +import numpy as np +import os +from openfermion import load_operator + +from tangelo.molecule_library import mol_H2_sto3g +from tangelo.toolboxes.operators.operators import QubitOperator +from tangelo.toolboxes.qubit_mappings import jordan_wigner, symmetry_conserving_bravyi_kitaev +from tangelo.toolboxes.ansatz_generator.vsqs import VSQS +from tangelo.linq import Simulator +from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit + +# For openfermion.load_operator function. +pwd_this_test = os.path.dirname(os.path.abspath(__file__)) + + +class VSQSTest(unittest.TestCase): + + def test_vsqs_set_var_params(self): + """Verify behavior of set_var_params for different inputs (list, numpy array). + MP2 have their own tests. + """ + + vsqs_ansatz = VSQS(mol_H2_sto3g) + + two_ones = np.ones((2,)) + + vsqs_ansatz.set_var_params([1., 1.]) + np.testing.assert_array_almost_equal(vsqs_ansatz.var_params, two_ones, decimal=6) + + vsqs_ansatz.set_var_params(np.array([1., 1.])) + np.testing.assert_array_almost_equal(vsqs_ansatz.var_params, two_ones, decimal=6) + + def test_vsqs_incorrect_number_var_params(self): + """Return an error if user provide incorrect number of variational + parameters. + """ + + vsqs_ansatz = VSQS(mol_H2_sto3g) + + self.assertRaises(ValueError, vsqs_ansatz.set_var_params, np.array([1., 1., 1., 1.])) + + def test_vsqs_H2(self): + """Verify closed-shell VSQS functionalities for H2.""" + + # Build circuit + vsqs_ansatz = VSQS(mol_H2_sto3g, intervals=3, time=3) + vsqs_ansatz.build_circuit() + + # Build qubit hamiltonian for energy evaluation + qubit_hamiltonian = jordan_wigner(mol_H2_sto3g.fermionic_hamiltonian) + + # Assert energy returned is as expected for given parameters + sim = Simulator() + vsqs_ansatz.update_var_params([0.66666667, 0.9698286, 0.21132472, 0.6465473]) + energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit) + self.assertAlmostEqual(energy, -1.1372701255155757, delta=1e-6) + + def test_vsqs_H4_doublecation(self): + """Verify closed-shell VSQS functionalities for H4 2+ by using saved qubit hamiltonian and initial hamiltonian""" + + var_params = [-2.53957674, 0.72683888, 1.08799500, 0.49836183, + -0.23020698, 0.93278630, 0.50591026, 0.50486903] + + # Build qubit hamiltonian for energy evaluation + qubit_hamiltonian = load_operator("mol_H4_doublecation_minao_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True) + initial_hamiltonian = load_operator("mol_H4_doublecation_minao_init_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True) + reference_state = get_reference_circuit(8, 2, "jw", up_then_down=True, spin=0) + + # Build circuit + vsqs_ansatz = VSQS(qubit_hamiltonian=qubit_hamiltonian, hini=initial_hamiltonian, reference_state=reference_state, + intervals=5, time=5, trotter_order=2) + vsqs_ansatz.build_circuit() + + # Assert energy returned is as expected for given parameters + sim = Simulator() + vsqs_ansatz.update_var_params(var_params) + energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit) + self.assertAlmostEqual(energy, -0.85425, delta=1e-4) + + def test_vsqs_H2_with_hnav(self): + """Verify closed-shell VSQS functionalities for H2 with navigator hamiltonian""" + navigator_hamiltonian = (QubitOperator('X0 Y1', 0.03632537110234512) + QubitOperator('Y0 X1', 0.03632537110234512) + + QubitOperator('Y0', 2.e-5) + QubitOperator('Y1', 2.e-5)) + + # Build qubit hamiltonian for energy evaluation + qubit_hamiltonian = symmetry_conserving_bravyi_kitaev(mol_H2_sto3g.fermionic_hamiltonian, 4, 2, False, 0) + + # Build circuit + vsqs_ansatz = VSQS(mol_H2_sto3g, intervals=2, time=1, mapping='scbk', up_then_down=True, trotter_order=2, + hnav=navigator_hamiltonian) + vsqs_ansatz.build_circuit() + + # Assert energy returned is as expected for given parameters + sim = Simulator() + vsqs_ansatz.update_var_params([0.50000001, -0.02494214, 3.15398767]) + energy = sim.get_expectation_value(qubit_hamiltonian, vsqs_ansatz.circuit) + self.assertAlmostEqual(energy, -1.1372701255155757, delta=1e-6) + + +if __name__ == "__main__": + unittest.main() diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py new file mode 100644 index 000000000..54bb492dd --- /dev/null +++ b/tangelo/toolboxes/ansatz_generator/vsqs.py @@ -0,0 +1,216 @@ +# Copyright 2021 Good Chemistry Company. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""This module defines the Variationally Scheduled Quantum Simulation class. It provides an +Adiabatic State Preparation (ASP) inspired ansatz as described in https://arxiv.org/abs/2003.09913.""" + +from copy import deepcopy + +import numpy as np +from openfermion import QubitOperator as ofQubitOperator + +from .ansatz import Ansatz +from .ansatz_utils import get_exponentiated_qubit_operator_circuit +from tangelo.toolboxes.operators import FermionOperator +from tangelo.linq import Circuit +from tangelo.toolboxes.qubit_mappings.mapping_transform import get_qubit_number, fermion_to_qubit_mapping +from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit + + +class VSQS(Ansatz): + """This class implements the Variationally Scheduled Quantum Simulator (VSQS) for state preparation as described in + https://arxiv.org/abs/2003.09913 + + Must supply either a molecule or a qubit_hamiltonian. If supplying a qubit_hamiltonian, must also supply + a reference_state Circuit and a hini QubitOperator. + + Args: + molecule (SecondQuantizedMolecule): The molecular system. Default: None + mapping (str): One of the support fermion to qubit mappings. Default : "JW" + up_then_down (bool): change basis ordering putting all spin up orbitals first, followed by all spin down. + Default: False (alternating spin up/down ordering) + intervals (int): The number of steps in the VSQS process. Must be greater than 1. Default: 2 + time (float): The propagation time. Default: 1 + qubit_hamiltonian (QubitOperator): The qubit hamiltonian to evolve. Default: None + reference_state (Circuit): The reference state for the propagation as defined by a Circuit. Mandatory if supplying + a qubit_hamiltonian. Default: None + hini (QubitOperator): The initial qubit Hamiltonian that one corresponds to the reference state. Mandatory if supplying + a qubit_hamiltonian. Default: None + hnav (QubitOperator): The navigator Hamiltonian. Default: None + """ + + def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2, time=1, qubit_hamiltonian=None, + reference_state=None, hini=None, hnav=None, trotter_order=1): + + self.up_then_down = up_then_down + self.mapping = mapping + if intervals > 1: + self.intervals = intervals + else: + raise ValueError("The number of intervals must be greater than 1.") + self.time = time + self.dt = self.time/self.intervals + self.reference_state = reference_state + self.trotter_order = trotter_order + + if molecule is None: + self.qubit_hamiltonian = qubit_hamiltonian + if not isinstance(hini, ofQubitOperator): + raise ValueError("When providing a qubit hamiltonian, an initial Hamiltonian must also be provided") + self.hini = hini + if not isinstance(reference_state, Circuit): + raise ValueError("Reference state must be provided when simulating a qubit hamiltonian directly") + self.reference_state = reference_state + else: + self.n_electrons = molecule.n_active_electrons + self.n_spinorbitals = int(molecule.n_sos) + self.n_qubits = get_qubit_number(mapping, self.n_spinorbitals) + self.spin = molecule.spin + self.qubit_hamiltonian = fermion_to_qubit_mapping(molecule.fermionic_hamiltonian, n_spinorbitals=self.n_spinorbitals, + n_electrons=self.n_electrons, mapping=self.mapping, + up_then_down=self.up_then_down, spin=self.spin) + self.hini = self.build_hini(molecule) if hini is None else hini + self.hfin = self.qubit_hamiltonian + self.hnav = hnav + + self.hfin_list, self.n_hfin = self.remove_constant_and_count_terms(self.qubit_hamiltonian) + self.hini_list, self.n_hini = self.remove_constant_and_count_terms(self.hini) + if self.hnav is None: + self.stride = 2 + self.n_hnav = 0 + else: + if isinstance(self.hnav, ofQubitOperator): + self.stride = 3 + self.hnav_list, self.n_hnav = self.remove_constant_and_count_terms(self.hnav) + else: + raise ValueError("Navigator Hamiltonian must be a QubitOperator") + + self.n_var_params = (intervals - 1) * self.stride + self.n_var_gates = (self.n_hini + self.n_hfin + self.n_hnav) * self.trotter_order + + self.var_params = None + self.circuit = None + + def remove_constant_and_count_terms(self, qu_op): + """count of non-zero terms in a QubitOperator and return the list of terms.items()""" + count = 0 + new_qu_op = deepcopy(qu_op) + for term in new_qu_op.terms.keys(): + if term: + count += 1 + else: + new_qu_op.terms[term] = 0 + new_qu_op.compress() + return list(new_qu_op.terms.items()), count + + def build_hini(self, molecule): + """Return the initial Hamiltonian (hini) composed of the one-body terms derived from the diagonal of Fock + matrix and one-body off-diagonal terms""" + fock = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_fock() @ molecule.mean_field.mo_coeff + h1 = molecule.mean_field.mo_coeff.T @ molecule.mean_field.get_hcore() @ molecule.mean_field.mo_coeff + hf_ferm = FermionOperator() + for i in range(self.n_spinorbitals//2): + for j in range(self.n_spinorbitals//2): + if i != j: + hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), h1[i, j]) + hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), h1[i, j]) + else: + hf_ferm += FermionOperator(((i * 2, 1), (j * 2, 0)), fock[i, j]) + hf_ferm += FermionOperator(((i * 2 + 1, 1), (j * 2 + 1, 0)), fock[i, j]) + return fermion_to_qubit_mapping(hf_ferm, mapping=self.mapping, n_spinorbitals=self.n_spinorbitals, n_electrons=self.n_electrons, + up_then_down=self.up_then_down, spin=self.spin) + + def set_var_params(self, var_params=None): + """Set values for the variational parameters. Default is linear interpolation.""" + if var_params is None: + var_params = self.init_params()[self.stride:self.n_var_params+self.stride] + + init_var_params = np.array(var_params) + if init_var_params.size == self.n_var_params: + self.var_params = var_params + else: + raise ValueError(f"Expected {self.n_var_params} variational parameters but received {init_var_params.size}.") + return var_params + + def update_var_params(self, var_params): + """Update the variational parameters in the circuit without rebuilding.""" + for i in range(self.intervals-1): + self._update_gate_params_for_qu_op(self.hini_list, self.n_var_gates * i, var_params[self.stride*i], self.n_hini) + self._update_gate_params_for_qu_op(self.hfin_list, self.n_var_gates * i + self.n_hini * self.trotter_order, + var_params[self.stride*i+1], self.n_hfin) + if self.hnav is not None: + self._update_gate_params_for_qu_op(self.hnav_list, self.n_var_gates * i + (self.n_hini + self.n_hfin) * self.trotter_order, + var_params[self.stride*i+2], self.n_hnav) + + def _update_gate_params_for_qu_op(self, qu_op, n_var_start, var_param, num_terms): + prefac = 2 / self.trotter_order * self.dt * var_param + for i, (_, coeff) in enumerate(qu_op): + self.circuit._variational_gates[n_var_start+i].parameter = prefac * coeff + if self.trotter_order == 2: + for i, (_, coeff) in enumerate(list(reversed(qu_op))): + self.circuit._variational_gates[n_var_start+num_terms+i].parameter = prefac * coeff + + def init_params(self): + """Generate initial parameters for the VSQS algorithm. + a_i = step*i, b_i=1-step*i, c_i= 1-step*i in_intervals/2""" + a = np.zeros(self.intervals+1) + b = np.zeros(self.intervals+1) + c = np.zeros(self.intervals+1) + a[0] = 1 + b[self.intervals] = 1 + step = 1/self.intervals + for i in range(1, self.intervals): + a[i] = (1 - i * step) + b[i] = (i * step) + all_params = np.zeros(self.stride * (self.intervals + 1)) + if self.hnav is None: + for i in range(self.intervals + 1): + all_params[2 * i] = a[i] + all_params[2 * i + 1] = b[i] + else: + c[0:self.intervals//2] = b[0:self.intervals//2] + c[self.intervals//2:self.intervals+1] = a[self.intervals//2:self.intervals+1] + for i in range(self.intervals + 1): + all_params[3 * i] = a[i] + all_params[3 * i + 1] = b[i] + all_params[3 * i + 2] = c[i] + return all_params + + def prepare_reference_state(self): + """Prepare a circuit generating the HF reference state.""" + return get_reference_circuit(n_spinorbitals=self.n_spinorbitals, n_electrons=self.n_electrons, mapping=self.mapping, + up_then_down=self.up_then_down, spin=self.spin) + + def build_circuit(self, var_params=None): + """Build the VSQS circuit by successive first-order trotterizations of hini, hfin and possibly hnav""" + reference_state_circuit = self.prepare_reference_state() if self.reference_state is None else self.reference_state + self.var_params = self.set_var_params(var_params) + + vsqs_circuit = get_exponentiated_qubit_operator_circuit(self.hini, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.hini_list) + for i in range(self.intervals-1): + vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hini, time=self.var_params[i * self.stride] * self.dt, variational=True, + trotter_order=self.trotter_order, pauli_order=self.hini_list) + vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hfin, time=self.var_params[i * self.stride + 1] * self.dt, variational=True, + trotter_order=self.trotter_order, pauli_order=self.hfin_list) + if self.hnav is not None: + vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hnav, time=self.var_params[i * self.stride + 2] * self.dt, variational=True, + trotter_order=self.trotter_order, pauli_order=self.hnav_list) + vsqs_circuit += get_exponentiated_qubit_operator_circuit(self.hfin, time=self.dt, trotter_order=self.trotter_order, pauli_order=self.hfin_list) + + if reference_state_circuit.size != 0: + self.circuit = reference_state_circuit + vsqs_circuit + else: + self.circuit = vsqs_circuit + + return self.circuit