diff --git a/tangelo/algorithms/projective/__init__.py b/tangelo/algorithms/projective/__init__.py new file mode 100644 index 000000000..a597fe080 --- /dev/null +++ b/tangelo/algorithms/projective/__init__.py @@ -0,0 +1,15 @@ +# 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. + +from .quantum_imaginary_time import QITESolver diff --git a/tangelo/algorithms/projective/quantum_imaginary_time.py b/tangelo/algorithms/projective/quantum_imaginary_time.py new file mode 100644 index 000000000..f6c47db83 --- /dev/null +++ b/tangelo/algorithms/projective/quantum_imaginary_time.py @@ -0,0 +1,315 @@ +# 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. + +"""Module that defines the Quantum Imaginary Time Algorithm (QITE) +""" +from copy import copy +import math + +from openfermion import FermionOperator as ofFermionOperator +import numpy as np + +from tangelo.toolboxes.ansatz_generator.ansatz_utils import trotterize +from tangelo.toolboxes.operators.operators import FermionOperator, QubitOperator +from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping +from tangelo.toolboxes.ansatz_generator._general_unitary_cc import uccgsd_generator as uccgsd_pool +from tangelo.toolboxes.operators import qubitop_to_qubitham +from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit +from tangelo.linq import Circuit, Simulator + + +class QITESolver: + """QITE class. This is an iterative algorithm that obtains a unitary operator + that approximates the imaginary time evolution of an initial state. + + Attributes: + molecule (SecondQuantizedMolecule): The molecular system. + dt (float): The imaginary time step size + min_de (float): Maximum energy change allowed before convergence. + max_cycles (int): Maximum number of iterations of QITE. + pool (func): Function that returns a list of FermionOperator. Each + element represents an excitation/operator that has an effect on the + total energy. + pool_args (tuple) : The arguments for the pool function given as a + tuple. + qubit_mapping (str): One of the supported qubit mapping identifiers. + qubit_hamiltonian (QubitOperator-like): Self-explanatory. + up_then_down (bool): Spin orbitals ordering. + n_spinorbitals (int): Self-explanatory. + n_electrons (int): Self-explanatory. + backend_options (dict): Backend options for the underlying QITE propagation + verbose (bool): Flag for verbosity of QITE. + """ + + def __init__(self, opt_dict): + + default_backend_options = {"target": None, "n_shots": None, "noise_model": None} + default_options = {"molecule": None, + "dt": 0.5, "max_cycles": 100, + "min_de": 1.e-7, + "pool": uccgsd_pool, + "pool_args": None, + "frozen_orbitals": "frozen_core", + "qubit_mapping": "jw", + "qubit_hamiltonian": None, + "up_then_down": False, + "n_spinorbitals": None, + "n_electrons": None, + "backend_options": default_backend_options, + "verbose": True} + + # Initialize with default values + self.__dict__ = default_options + # Overwrite default values with user-provided ones, if they correspond to a valid keyword + for k, v in opt_dict.items(): + if k in default_options: + setattr(self, k, v) + else: + raise KeyError(f"Keyword :: {k}, not available in {self.__class__.__name__}") + + # Raise error/warnings if input is not as expected. Only a single input + # must be provided to avoid conflicts. + if not (bool(self.molecule) ^ bool(self.qubit_hamiltonian)): + raise ValueError("A molecule OR qubit Hamiltonian object must be provided when instantiating " + f"{self.__class__.__name__}.") + + if self.qubit_hamiltonian: + if not (self.n_spinorbitals and self.n_electrons): + raise ValueError("Expecting the number of spin-orbitals (n_spinorbitals) and the number of " + "electrons (n_electrons) with a qubit_hamiltonian.") + + self.iteration = 0 + self.energies = list() + + self.circuit_list = list() + self.final_energy = None + self.final_circuit = None + self.final_statevector = None + + self.backend = None + + def prepare_reference_state(self): + """Returns circuit preparing the reference state of the ansatz (e.g + prepare reference wavefunction with HF, multi-reference state, etc). + These preparations must be consistent with the transform used to obtain + the qubit operator. + """ + + return get_reference_circuit(n_spinorbitals=self.n_spinorbitals, + n_electrons=self.n_electrons, + mapping=self.qubit_mapping, + up_then_down=self.up_then_down, + spin=0) + + def build(self): + """Builds the underlying objects required to run the QITE algorithm.""" + + # Building molecule data with a pyscf molecule. + if self.molecule: + + self.n_spinorbitals = self.molecule.n_active_sos + self.n_electrons = self.molecule.n_active_electrons + + # Compute qubit hamiltonian for the input molecular system + qubit_op = fermion_to_qubit_mapping(fermion_operator=self.molecule.fermionic_hamiltonian, + mapping=self.qubit_mapping, + n_spinorbitals=self.n_spinorbitals, + n_electrons=self.n_electrons, + up_then_down=self.up_then_down, + spin=0) + + self.qubit_hamiltonian = qubitop_to_qubitham(qubit_op, self.qubit_mapping, self.up_then_down) + + # Getting the pool of operators for the ansatz. If more functionalities + # are added, this part must be modified and generalized. + if self.pool_args is None: + if self.pool == uccgsd_pool: + self.pool_args = (self.n_spinorbitals,) + else: + raise KeyError("pool_args must be defined if using own pool function") + # Check if pool function returns a QubitOperator or FermionOperator and populate variables + pool_list = self.pool(*self.pool_args) + if isinstance(pool_list[0], QubitOperator): + self.pool_type = "qubit" + self.full_pool_operators = pool_list + elif isinstance(pool_list[0], (FermionOperator, ofFermionOperator)): + self.pool_type = "fermion" + self.fermionic_operators = pool_list + self.full_pool_operators = [fermion_to_qubit_mapping(fermion_operator=fi, + mapping=self.qubit_mapping, + n_spinorbitals=self.n_spinorbitals, + n_electrons=self.n_electrons, + up_then_down=self.up_then_down) for fi in self.fermionic_operators] + else: + raise ValueError("pool function must return either QubitOperator or FermionOperator") + + # Cast all coefs to floats (rotations angles are real). + for qubit_op in self.full_pool_operators: + for term, coeff in qubit_op.terms.items(): + qubit_op.terms[term] = math.copysign(1., coeff.imag) + + # Remove duplicates and only select terms with odd number of Y gates for all mappings except JKMN + if self.qubit_mapping.upper() != "JKMN": + reduced_pool_terms = set() + for qubit_op in self.full_pool_operators: + for term in qubit_op.terms: + count_y = str(term).count("Y") + if count_y % 2 == 1: + reduced_pool_terms.add(term) + else: + reduced_pool_terms = set() + for qubit_op in self.full_pool_operators: + for term in qubit_op.terms.keys(): + if term: + reduced_pool_terms.add(term) + + # Generated list of pool_operators and full pool operator. + self.pool_operators = [QubitOperator(term) for term in reduced_pool_terms] + self.pool_qubit_op = QubitOperator() + for term in self.pool_operators: + self.pool_qubit_op += term + + self.qubit_operator = self.qubit_hamiltonian.to_qubitoperator() + + # Obtain all qubit terms that need to be measured + self.pool_h = [element*self.qubit_operator for element in self.pool_operators] + self.pool_pool = [[element1*element2 for element2 in self.pool_operators] for element1 in self.pool_operators] + + # Obtain initial state preparation circuit + self.circuit_list.append(self.prepare_reference_state()) + self.final_circuit = copy(self.circuit_list[0]) + + # Quantum circuit simulation backend options + self.backend = Simulator(target=self.backend_options["target"], n_shots=self.backend_options["n_shots"], + noise_model=self.backend_options["noise_model"]) + + self.use_statevector = self.backend.statevector_available and self.backend._noise_model is None + + def simulate(self): + """Performs the QITE cycles. Each iteration, a linear system is + solved to obtain the next unitary. The system to be solved can be found in + section 3.5 of https://arxiv.org/pdf/2108.04413.pdf + + Returns: + float: final energy after obtaining running QITE + """ + + # Construction of the circuit. self.max_cycles terms are added, unless + # the energy change is less than self.min_de. + if self.use_statevector: + self.update_statevector(self.backend, self.circuit_list[0]) + self.final_energy = self.energy_expectation(self.backend) + self.energies.append(self.final_energy) + while self.iteration < self.max_cycles: + self.iteration += 1 + if self.verbose: + print(f"Iteration {self.iteration} of QITE with starting energy {self.final_energy}") + + suv, bu = self.calculate_matrices(self.backend, self.final_energy) + + alphas_array = np.linalg.solve(suv.real, bu.real) + # convert to dictionary with key as first (only) term of each pool_operator and value self.dt * alphas_array[i] + alphas_dict = {next(iter(qu_op.terms)): self.dt * alphas_array[i] for i, qu_op in enumerate(self.pool_operators)} + next_circuit = trotterize(self.pool_qubit_op, alphas_dict, trotter_order=1, n_trotter_steps=1) + + self.circuit_list.append(next_circuit) + self.final_circuit += next_circuit + + if self.use_statevector: + self.update_statevector(self.backend, self.circuit_list[self.iteration]) + + new_energy = self.energy_expectation(self.backend) + self.energies.append(new_energy) + + if abs(new_energy - self.final_energy) < self.min_de and self.iteration > 1: + self.final_energy = new_energy + break + self.final_energy = new_energy + + if self.verbose: + print(f"Final energy of QITE is {self.final_energy}") + + return self.energies[-1] + + def update_statevector(self, backend: Simulator, next_circuit: Circuit): + r"""Update self.final_statevector by propagating with next_circuit using backend + + Args: + Simulator: the backend to use for the statevector update + Circuit: The circuit to apply to the statevector + """ + _, self.final_statevector = backend.simulate(next_circuit, + return_statevector=True, + initial_statevector=self.final_statevector) + + def calculate_matrices(self, backend: Simulator, new_energy: float): + r"""Calculated matrix elements for imaginary time evolution. + The matrices are defined in section 3.5 of https://arxiv.org/pdf/2108.04413.pdf + + Args: + backend (Simulator): the backend from which the matrices are generated + new_energy (float): the current energy_expectation of the Hamiltonian + + Returns: + matrix float: The expectation values <\psi| pu^+ pv |\psi> + array float: The expecation values <\psi| pu^+ H |\psi> + """ + + circuit = Circuit(n_qubits=self.final_circuit.width) if self.use_statevector else self.final_circuit + + ndeltab = np.sqrt(1 - 2 * self.dt * new_energy) + prefac = -1j/ndeltab + bu = [prefac*backend.get_expectation_value(element, circuit, initial_statevector=self.final_statevector) + for element in self.pool_h] + bu = np.array(bu) + pool_size = len(self.pool_h) + suv = np.zeros((pool_size, pool_size), dtype=complex) + for u in range(pool_size): + for v in range(u+1, pool_size): + suv[u, v] = backend.get_expectation_value(self.pool_pool[u][v], + circuit, + initial_statevector=self.final_statevector) + suv[v, u] = suv[u, v] + + return suv, bu + + def energy_expectation(self, backend: Simulator): + """Estimate energy using the self.final_circuit, qubit hamiltonian and compute + backend. + + Args: + backend (Simulator): the backend one computes the energy expectation with + + Returns: + float: energy computed by the backend + """ + circuit = Circuit(n_qubits=self.final_circuit.width) if self.use_statevector else self.final_circuit + energy = backend.get_expectation_value(self.qubit_hamiltonian.to_qubitoperator(), + circuit, + initial_statevector=self.final_statevector) + return energy + + def get_resources(self): + """Returns resources currently used in underlying state preparation i.e. self.final_circuit + the number of pool operators, and the size of qubit_hamiltonian + + Returns: + dict: Dictionary of various quantum resources required""" + resources = dict() + resources["qubit_hamiltonian_terms"] = len(self.qubit_hamiltonian.terms) + resources["pool_size"] = len(self.pool_operators) + resources["circuit_width"] = self.final_circuit.width + resources["circuit_gates"] = self.final_circuit.size + resources["circuit_2qubit_gates"] = self.final_circuit.counts.get("CNOT", 0) + return resources diff --git a/tangelo/algorithms/projective/tests/__init__.py b/tangelo/algorithms/projective/tests/__init__.py new file mode 100644 index 000000000..532746351 --- /dev/null +++ b/tangelo/algorithms/projective/tests/__init__.py @@ -0,0 +1,13 @@ +# 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. diff --git a/tangelo/algorithms/projective/tests/test_qite.py b/tangelo/algorithms/projective/tests/test_qite.py new file mode 100644 index 000000000..c9f909bbc --- /dev/null +++ b/tangelo/algorithms/projective/tests/test_qite.py @@ -0,0 +1,158 @@ +# Copyright 2021 Good Chemsitry 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 + +from tangelo.algorithms.projective.quantum_imaginary_time import QITESolver +from tangelo.molecule_library import mol_H2_sto3g, mol_H4_sto3g +from tangelo.linq.noisy_simulation import NoiseModel + + +class QITESolverTest(unittest.TestCase): + + def test_instantiation_qite(self): + """Try instantiating QITESolver with basic input.""" + + options = {"molecule": mol_H2_sto3g, "qubit_mapping": "jw"} + QITESolver(options) + + def test_instantiation_qite_incorrect_keyword(self): + """Instantiating with an incorrect keyword should return an error """ + + options = {"molecule": mol_H2_sto3g, "qubit_mapping": "jw", "dummy": True} + self.assertRaises(KeyError, QITESolver, options) + + def test_instantiation_qite_missing_molecule(self): + """Instantiating with no molecule should return an error.""" + + options = {"qubit_mapping": "jw"} + self.assertRaises(ValueError, QITESolver, options) + + def test_simulate_h2_noisy(self): + """Run QITE on H2 molecule with bk qubit mapping and an empty noise model for 1 cycle. + Result should be lower than mean field energy. + """ + + backend_options = {"target": None, "n_shots": 10000, "noise_model": NoiseModel()} + + qite_options = {"molecule": mol_H2_sto3g, "qubit_mapping": "scbk", + "verbose": True, "backend_options": backend_options, + "max_cycles": 1, "up_then_down": True} + qite_solver = QITESolver(qite_options) + qite_solver.build() + + energy = qite_solver.simulate() + self.assertTrue(energy < mol_H2_sto3g.mf_energy) + + def test_simulate_h2(self): + """Run QITE on H2 molecule, with JW qubit mapping and exact simulator + """ + + qite_options = {"molecule": mol_H2_sto3g, "qubit_mapping": "jw", + "verbose": False, "up_then_down": True} + qite_solver = QITESolver(qite_options) + qite_solver.build() + + energy = qite_solver.simulate() + self.assertAlmostEqual(energy, -1.137270422018, delta=1e-4) + + def test_resources_h2(self): + """Test get_resources funtion for QITE on H2 molecule, with JW qubit mapping. + """ + + qite_options = {"molecule": mol_H2_sto3g, "qubit_mapping": "jw", + "verbose": False} + qite_solver = QITESolver(qite_options) + qite_solver.build() + resources = qite_solver.get_resources() + self.assertEqual(resources["qubit_hamiltonian_terms"], 15) + self.assertEqual(resources["pool_size"], 20) + + def test_mapping_BK(self): + """Test that BK mapping recovers the expected result for the example of H2. + """ + qite_options = {"molecule": mol_H2_sto3g, "verbose": False, + "qubit_mapping": "bk"} + + qite_solver = QITESolver(qite_options) + qite_solver.build() + energy = qite_solver.simulate() + + energy_target = -1.137270 + self.assertAlmostEqual(energy, energy_target, places=5) + + def test_mapping_JKMN(self): + """Test that JKMN mapping recovers the expected result for the example of H2. + """ + qite_options = {"molecule": mol_H2_sto3g, "verbose": False, + "qubit_mapping": "JKMN"} + + qite_solver = QITESolver(qite_options) + qite_solver.build() + energy = qite_solver.simulate() + + energy_target = -1.137270 + self.assertAlmostEqual(energy, energy_target, places=5) + + def test_mapping_scBK(self): + """Test that scBK mapping recovers the expected result for the example of H2. + """ + qite_options = {"molecule": mol_H2_sto3g, "verbose": False, + "qubit_mapping": "scbk", "up_then_down": True} + + qite_solver = QITESolver(qite_options) + qite_solver.build() + energy = qite_solver.simulate() + + energy_target = -1.137270 + self.assertAlmostEqual(energy, energy_target, places=5) + + def test_spin_reorder_equivalence(self): + """Test that re-ordered spin input (all up followed by all down) returns + the same optimized energy result for both JW and BK mappings. + """ + qite_options = {"molecule": mol_H2_sto3g, "up_then_down": True, + "verbose": False, "qubit_mapping": "jw"} + + qite_solver_jw = QITESolver(qite_options) + qite_solver_jw.build() + energy_jw = qite_solver_jw.simulate() + + qite_options["qubit_mapping"] = "bk" + qite_solver_bk = QITESolver(qite_options) + qite_solver_bk.build() + energy_bk = qite_solver_bk.simulate() + + energy_target = -1.137270 + self.assertAlmostEqual(energy_jw, energy_target, places=5) + self.assertAlmostEqual(energy_bk, energy_target, places=5) + + def test_simulate_h4_frozen_orbitals(self): + """Run QITE on H4 molecule, with UCCSD ansatz, JW qubit mapping, initial + parameters, exact simulator. First (occupied) and last (virtual) + orbitals are frozen. + """ + mol_H4_sto3g_frozen = mol_H4_sto3g.freeze_mos([0, 3], inplace=False) + + qite_options = {"molecule": mol_H4_sto3g_frozen, "qubit_mapping": "jw", + "verbose": False} + qite_solver = QITESolver(qite_options) + qite_solver.build() + + energy = qite_solver.simulate() + self.assertAlmostEqual(energy, -1.8943598012229799, delta=1e-5) + + +if __name__ == "__main__": + unittest.main() diff --git a/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py b/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py index c4531a5d6..94a3c5d63 100644 --- a/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py +++ b/tangelo/toolboxes/ansatz_generator/adapt_ansatz.py @@ -18,7 +18,7 @@ from tangelo.linq import Circuit from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit -from tangelo.toolboxes.ansatz_generator.ansatz_utils import pauliword_to_circuit +from tangelo.toolboxes.ansatz_generator.ansatz_utils import exp_pauliword_to_gates from tangelo.toolboxes.ansatz_generator.ansatz import Ansatz from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping @@ -118,7 +118,7 @@ def build_circuit(self, var_params=None): for op in self.operators: for pauli_term in op.get_operators(): pauli_tuple = list(pauli_term.terms.keys())[0] - adapt_circuit += pauliword_to_circuit(pauli_tuple, 0.1) + adapt_circuit += exp_pauliword_to_gates(pauli_tuple, 0.1) adapt_circuit = Circuit(adapt_circuit) if adapt_circuit.size != 0: @@ -156,6 +156,6 @@ def add_operator(self, pauli_operator, ferm_operator=None): self._var_params_prefactor += [math.copysign(1., coeff)] pauli_tuple = list(pauli_term.terms.keys())[0] - new_operator = Circuit(pauliword_to_circuit(pauli_tuple, 0.1)) + new_operator = Circuit(exp_pauliword_to_gates(pauli_tuple, 0.1)) self.circuit += new_operator diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py index b2b546fda..a9a2a4c16 100644 --- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py +++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py @@ -17,8 +17,17 @@ facilitate the assembly of ansatz quantum circuits. """ +from copy import deepcopy +from itertools import combinations + import numpy as np +from openfermion.ops import FermionOperator as ofFermionOperator +from openfermion.ops import InteractionOperator as ofInteractionOperator +from openfermion.ops import QubitOperator as ofQubitOperator + from tangelo.linq import Circuit, Gate +from tangelo.toolboxes.operators import FermionOperator, QubitOperator +from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping, get_fermion_operator def pauli_op_to_gate(index, op, inverse=False): @@ -32,9 +41,18 @@ def pauli_op_to_gate(index, op, inverse=False): return Gate("RX", index, parameter=angle) if not inverse else Gate("RX", index, parameter=-angle+4*np.pi) -def pauliword_to_circuit(pauli_word, coef, variational=True): - """Generates a quantum circuit corresponding to the pauli word, as described - in Whitfield 2010 (https://arxiv.org/pdf/1001.3855.pdf). +def exp_pauliword_to_gates(pauli_word, coef, variational=True, control=None): + """Generate a list of Gate objects corresponding to the exponential of a pauli word. + The process is described in Whitfield 2010 https://arxiv.org/pdf/1001.3855.pdf + + Args: + pauli_word (tuple): Openfermion-like tuple that generates a pauli_word to exponentiate + coef (float): The coefficient in the exponentiation + variational (bool): When creating the Gate objects, label the (controlled-)Rz gate as variational + control (integer): The control qubit label + + Returns: + list: list of Gate objects that represents the exponentiation of the pauli word. """ gates = [] @@ -49,7 +67,10 @@ def pauliword_to_circuit(pauli_word, coef, variational=True): gates += cnot_ladder_gates angle = 2.*coef if coef >= 0. else 4*np.pi+2*coef - gates += [Gate("RZ", target=indices[-1], parameter=angle, is_variational=variational)] + if control is None: + gates += [Gate("RZ", target=indices[-1], parameter=angle, is_variational=variational)] + else: + gates += [Gate("CRZ", target=indices[-1], control=control, parameter=angle)] gates += cnot_ladder_gates[::-1] @@ -59,3 +80,317 @@ def pauliword_to_circuit(pauli_word, coef, variational=True): gates += [pauli_op_to_gate(index, op, inverse=True)] return gates + + +def get_exponentiated_qubit_operator_circuit(qubit_op, time=1., variational=False, trotter_order=1, control=None, + return_phase=False): + """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 + + Args: + qubit_op (QubitOperator): qubit hamiltonian to exponentiate + time (float or dict): The time to evolve the whole system or individiual times for each + term in the operator. If a dictionary, must have keys that have a matching key in qubit_op.terms + 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 + + 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 trotter_order > 2: + raise ValueError(f"Trotter order of >2 is not supported currently in Tangelo.") + prefactor = 1/2 if trotter_order == 2 else 1 + + if isinstance(time, (float, np.floating, np.integer, int)): + evolve_time = {term: prefactor*time for term in qubit_op.terms.keys()} + elif isinstance(time, dict): + if time.keys() == qubit_op.terms.keys(): + evolve_time = {term: prefactor*etime for term, etime in time.items()} + else: + raise ValueError(f"The keys in the time dictionary do not match the keys in qubit_op.terms") + else: + raise ValueError(f"time must be a float or a dictionary") + + phase = 1. + exp_pauli_word_gates = list() + for i in range(trotter_order): + if i == 0: + pauli_words.reverse() + for pauli_word, coef in pauli_words: + if pauli_word: # identity terms do not contribute to evolution outside of a phase + if abs(np.real(coef)*evolve_time[pauli_word]) > 1.e-10: + exp_pauli_word_gates += exp_pauliword_to_gates(pauli_word, + np.real(coef)*evolve_time[pauli_word], + variational=variational, + control=control) + else: + if control is None: + phase *= np.exp(-1j * coef * evolve_time[pauli_word]) + else: + exp_pauli_word_gates += [Gate("PHASE", target=control, parameter=-np.real(coef)*evolve_time[pauli_word])] + + return_value = (Circuit(exp_pauli_word_gates), phase) if return_phase else Circuit(exp_pauli_word_gates) + return return_value + + +def trotterize(operator, time=1., n_trotter_steps=1, trotter_order=1, variational=False, + mapping_options=dict(), control=None, return_phase=False): + """Generate the circuit that represents time evolution of an operator. + This circuit is generated as a trotterization of a qubit operator which is either the input + or mapped from the given fermion operator. + + Args: + operator (QubitOperator or FermionOperator): operator to time evolve + time (float or dict): The time to evolve the whole system or individiual times for each + term in the operator. If a dict, each key must match the keys in operator.terms + variational (bool): whether the coefficients are variational + trotter_order (int): order of trotter approximation, 1 or 2 supported + n_trotter_steps (int): The number of different time steps taken for total time t + mapping_options (dict): Defines the desired Fermion->Qubit mapping + Default values:{"up_then_down": False, "qubit_mapping": "jw", "n_spinorbitals": None, + "n_electrons": None} + control (int): The label for the control Qubit of the time-evolution + return_phase (bool): If return_phase is True, the global phase of the time-evolution will be returned + + Returns: + Circuit: circuit corresponding to time evolution of the operator + float: the global phase not included in the circuit if return_phase=True else not included + """ + + if isinstance(operator, ofFermionOperator): + options = {"up_then_down": False, "qubit_mapping": "jw", "n_spinorbitals": None, "n_electrons": None} + # Overwrite default values with user-provided ones, if they correspond to a valid keyword + for k, v in mapping_options.items(): + if k in options: + options[k] = v + else: + raise KeyError(f"Keyword :: {k}, not a valid fermion to qubit mapping option") + if isinstance(time, (float, np.floating, int, np.integer)): + evolve_time = {term: time for term in operator.terms.keys()} + elif isinstance(time, dict): + if time.keys() == operator.terms.keys(): + evolve_time = deepcopy(time) + else: + raise ValueError(f"keys of time do not match keys of operator.terms") + else: + raise ValueError("time must be a float or dictionary of floats") + new_operator = FermionOperator() + for term in operator.terms: + new_operator += FermionOperator(term, operator.terms[term]*evolve_time[term]/n_trotter_steps) + qubit_op = fermion_to_qubit_mapping(fermion_operator=new_operator, + mapping=options["qubit_mapping"], + n_spinorbitals=options["n_spinorbitals"], + n_electrons=options["n_electrons"], + up_then_down=options["up_then_down"]) + circuit, phase = get_exponentiated_qubit_operator_circuit(qubit_op, + time=1., # time is already included + trotter_order=trotter_order, + variational=variational, + control=control, + return_phase=True) + + elif isinstance(operator, (ofQubitOperator)): + qubit_op = deepcopy(operator) + if isinstance(time, float): + evolve_time = time / n_trotter_steps + elif isinstance(time, dict): + if time.keys() == operator.terms.keys(): + evolve_time = {term: etime / n_trotter_steps for term, etime in time.items()} + else: + raise ValueError(f"time dictionary and operator.terms dictionary have different keys.") + else: + raise ValueError(f"time must be a float or a dictionary of floats") + circuit, phase = get_exponentiated_qubit_operator_circuit(qubit_op, + time=evolve_time, + trotter_order=trotter_order, + variational=variational, + control=control, + return_phase=True) + else: + raise ValueError("Only FermionOperator or QubitOperator allowed") + + if n_trotter_steps == 1: + return_value = (circuit, phase) if return_phase else circuit + else: + final_circuit = deepcopy(circuit) + final_phase = deepcopy(phase) + for i in range(1, n_trotter_steps): + final_circuit += circuit + final_phase *= phase + return_value = (final_circuit, final_phase) if return_phase else circuit + return return_value + + +def append_qft_rotations_gates(gate_list, qubit_list, prefac=1): + """Appends the list of gates required for a quantum fourier transform to a gate list. + + Args: + gate_list (list): List of Gate elements + qubit_list (list): List of integers for which the qft operations are performed + + Returns: + list: List of Gate objects for rotation portion of qft circuit appended to gate_list""" + n = len(qubit_list) + if n == 0: + return gate_list + n -= 1 + gate_list += [Gate("H", target=qubit_list[n])] + for i, qubit in enumerate(qubit_list[:n]): + gate_list += [Gate("CPHASE", control=qubit, target=qubit_list[n], parameter=prefac*np.pi/2**(n-i))] + + append_qft_rotations_gates(gate_list, qubit_list[:n], prefac=prefac) + + +def swap_registers(gate_list, qubit_list): + """Function to swap register order. + Args: + gate_list (list): List of Gate + qubit_list (list): List of integers for the locations of the qubits + + Result: + list: The Gate operations that swap the register order""" + n = len(qubit_list) + for qubit_index in range(n//2): + gate_list += [Gate("SWAP", target=[qubit_list[qubit_index], qubit_list[n - qubit_index - 1]])] + return gate_list + + +def get_qft_circuit(qubits, n_qubits=None, inverse=False, swap=True): + """Returns the QFT or iQFT circuit given a list of qubits to act on. + + Args: + qubits (int or list): The list of qubits to apply the QFT circuit to. If an integer, + the operation is applied to the [0,...,qubits-1] qubits + n_qubits: Argument to initialize a Circuit with the desired number of qubits. + inverse (bool): If True, the inverse QFT is applied. If False, QFT is applied + swap (bool): Whether to apply swap to the registers. + + Returns: + Circuit: The circuit that applies QFT or iQFT to qubits + """ + + if isinstance(qubits, int): + qubit_list = list(range(qubits)) + elif isinstance(qubits, list): + qubit_list = qubits + else: + raise KeyError("qubits must be an int or list of ints") + + swap_gates = list() + if swap: + swap_registers(swap_gates, qubit_list) + + qft_gates = list() + if inverse: + append_qft_rotations_gates(qft_gates, qubit_list, prefac=-1) + qft_gates = [gate for gate in reversed(qft_gates)] + qft_gates = swap_gates + qft_gates + else: + append_qft_rotations_gates(qft_gates, qubit_list) + qft_gates += swap_gates + + return Circuit(qft_gates, n_qubits=n_qubits) + + +def controlled_pauliwords(qubit_op, control, n_qubits=None): + """Takes a qubit operator and returns controlled-pauliword circuits for each term as a list. + + Args: + qubit_op (QubitOperator): The qubit operator with pauliwords to generate circuits for + control (int): The index of the control qubit + n_qubits (int): When generating each Circuit, create with n_qubits size + + Returns: + list: List of controlled-pauliword Circuit for each pauliword in the qubit_op + """ + pauli_words = qubit_op.terms.items() + + pauliword_circuits = list() + for (pauli_word, _) in pauli_words: + gates = [Gate(name="C"+op, target=index, control=control) for index, op in pauli_word] + pauliword_circuits.append(Circuit(gates, n_qubits=n_qubits)) + return pauliword_circuits + + +def controlled_swap_to_XX_gates(c, n1, n2): + """Equivalent decomposition of controlled swap into 1-qubit gates and XX 2-qubit gate. + + This is useful for IonQ experiments as the native two-qubit gate is the XX Ising coupling. + + Args: + c (int): control qubit + n1 (int): first target qubit + n2 (int): second target qubit + + Returns: + list: List of Gate that applies controlled swap operation + """ + gates = [Gate("RY", target=c, parameter=7*np.pi/2.), + Gate("RZ", target=n1, parameter=7*np.pi/2.), + Gate("XX", target=[n1, n2], parameter=5*np.pi/2.), + Gate("RZ", target=n1, parameter=7*np.pi/4.), + Gate("RZ", target=n2, parameter=3*np.pi/4.), + Gate("RY", target=n1, parameter=np.pi/2.), + Gate("XX", target=[c, n2], parameter=7*np.pi/2.), + Gate("RY", target=n2, parameter=11*np.pi/4), + Gate("XX", target=[n1, n2], parameter=7*np.pi/2.), + Gate("XX", target=[c, n1], parameter=np.pi/4.), + Gate("RZ", target=n2, parameter=np.pi/4), + Gate("XX", target=[c, n2], parameter=5*np.pi/2), + Gate("RY", target=c, parameter=5*np.pi/2), + Gate("RZ", target=n1, parameter=5*np.pi/2), + Gate("RY", target=n2, parameter=7*np.pi/4), + Gate("XX", target=[n1, n2], parameter=7*np.pi/2), + Gate("RY", target=n1, parameter=np.pi/2), + Gate("RZ", target=c, parameter=11*np.pi/4)] + return gates + + +def derangement_circuit(qubit_list, control=None, n_qubits=None, decomp=None): + """Returns the (controlled-)derangement circuit for multiple copies of a state + + Args: + qubit_list (list of list(int)): Each item in the list is a list of qubit registers for each copy. The length of + each list of qubit registers must be the same. + For example [[1, 2], [3, 4]] applies controlled-swaps between equivalent states located on qubits [1, 2] and [3, 4] + control (int): The control register to be measured. + n_qubits (int): The number of qubits in the circuit. + decomp (str): Use the decomposed controlled-swap into 1-qubit gates and a certain 2-qubit gate listed below. + "XX": 2-qubit gate is XX + + Returns: + Circuit: The derangement circuit + """ + if decomp is not None and decomp not in ["XX"]: + raise ValueError(f"{decomp} is not a valid controlled swap decomposition") + + num_copies = len(qubit_list) + if num_copies == 1: + return Circuit(n_qubits=n_qubits) + else: + rho_range = len(qubit_list[0]) + for i in range(1, num_copies): + if len(qubit_list[i]) != rho_range: + raise ValueError("All copies must have the same number of qubits") + gate_list = list() + if control is None: + for copy1, copy2 in combinations(range(num_copies), 2): + for rhoi in range(rho_range): + gate_list += [Gate("SWAP", target=[qubit_list[copy1][rhoi], qubit_list[copy2][rhoi]])] + else: + for copy1, copy2 in combinations(range(num_copies), 2): + for rhoi in range(rho_range): + if decomp == "XX": + gate_list += controlled_swap_to_XX_gates(control, + qubit_list[copy1][rhoi], + qubit_list[copy2][rhoi]) + else: + gate_list += [Gate("CSWAP", + target=[qubit_list[copy1][rhoi], qubit_list[copy2][rhoi]], + control=control)] + + return Circuit(gate_list, n_qubits=n_qubits) diff --git a/tangelo/toolboxes/ansatz_generator/qcc.py b/tangelo/toolboxes/ansatz_generator/qcc.py index c6be5852e..737da5975 100755 --- a/tangelo/toolboxes/ansatz_generator/qcc.py +++ b/tangelo/toolboxes/ansatz_generator/qcc.py @@ -41,7 +41,7 @@ fermion_to_qubit_mapping from tangelo.linq import Circuit from .ansatz import Ansatz -from .ansatz_utils import pauliword_to_circuit +from .ansatz_utils import exp_pauliword_to_gates from ._qubit_mf import init_qmf_from_hf, get_qmf_circuit, purify_qmf_state from ._qubit_cc import construct_dis @@ -208,7 +208,7 @@ def build_circuit(self, var_params=None): pauli_words = sorted(qubit_op.terms.items(), key=lambda x: len(x[0])) pauli_words_gates = [] for i, (pauli_word, coef) in enumerate(pauli_words): - pauli_words_gates += pauliword_to_circuit(pauli_word, coef) + pauli_words_gates += exp_pauliword_to_gates(pauli_word, coef) self.pauli_to_angles_mapping[pauli_word] = i self.qcc_circuit = Circuit(pauli_words_gates) self.circuit = self.qmf_circuit + self.qcc_circuit if self.qmf_circuit.size != 0\ diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py b/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py new file mode 100644 index 000000000..10068e235 --- /dev/null +++ b/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py @@ -0,0 +1,301 @@ +import unittest + +from scipy.linalg import expm +import numpy as np +from numpy.linalg import eigh +from openfermion import get_sparse_operator + +from tangelo.linq import Simulator, Circuit, Gate +from tangelo.toolboxes.operators import FermionOperator, QubitOperator +from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping +from tangelo.molecule_library import mol_H4_sto3g +from tangelo.linq.tests.test_simulator import assert_freq_dict_almost_equal +from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit +from tangelo.toolboxes.ansatz_generator.ansatz_utils import trotterize, get_qft_circuit +from tangelo.toolboxes.ansatz_generator.ansatz_utils import controlled_swap_to_XX_gates +from tangelo.toolboxes.ansatz_generator.ansatz_utils import derangement_circuit, controlled_pauliwords + +# Initiate simulators, Use cirq as it has the same ordering for statevectors as openfermion does for Hamiltonians +# This is important when converting the openfermion QubitOperator toarray(), propagating exactly and comparing +# to the statevector output of the simulator. All other simulators will produce the same statevector values but +# in a different order (i.e. msq_first instead of lsq_first) +sim = Simulator(target="cirq") + +fermion_operator = mol_H4_sto3g._get_fermionic_hamiltonian() + + +class ansatz_utils_Test(unittest.TestCase): + + def test_trotterize_fermion_input(self): + """ Verify that the time evolution is correct for different mappings and a fermionic + hamiltonian input + """ + + time = 0.2 + for mapping in ["jw", "bk", "scbk"]: + reference_circuit = get_reference_circuit(n_spinorbitals=mol_H4_sto3g.n_active_sos, + n_electrons=mol_H4_sto3g.n_active_electrons, + mapping=mapping, + up_then_down=True) + _, refwave = sim.simulate(reference_circuit, return_statevector=True) + + qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=fermion_operator, + mapping=mapping, + n_spinorbitals=mol_H4_sto3g.n_active_sos, + n_electrons=mol_H4_sto3g.n_active_electrons, + up_then_down=True) + + ham_mat = get_sparse_operator(qubit_hamiltonian).toarray() + evolve_exact = expm(-1j * time * ham_mat) @ refwave + + options = {"up_then_down": True, + "qubit_mapping": mapping, + "n_spinorbitals": mol_H4_sto3g.n_active_sos, + "n_electrons": mol_H4_sto3g.n_active_electrons} + tcircuit, phase = trotterize(fermion_operator, trotter_order=1, n_trotter_steps=1, time=time, + mapping_options=options, return_phase=True) + _, wavefunc = sim.simulate(tcircuit, return_statevector=True, initial_statevector=refwave) + wavefunc *= phase + overlap = np.dot(np.conj(evolve_exact), wavefunc) + self.assertAlmostEqual(overlap, 1.0, delta=1e-3) + + def test_trotterize_qubit_input(self): + """ Verify that the time evolution is correct for different mappings and a qubit_hamiltonian input""" + + time = 0.2 + for mapping in ["jw", "bk", "scbk"]: + reference_circuit = get_reference_circuit(n_spinorbitals=mol_H4_sto3g.n_active_sos, + n_electrons=mol_H4_sto3g.n_active_electrons, + mapping=mapping, + up_then_down=True) + _, refwave = sim.simulate(reference_circuit, return_statevector=True) + + qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=fermion_operator, + mapping=mapping, + n_spinorbitals=mol_H4_sto3g.n_active_sos, + n_electrons=mol_H4_sto3g.n_active_electrons, + up_then_down=True) + + ham_mat = get_sparse_operator(qubit_hamiltonian).toarray() + evolve_exact = expm(-1j * time * ham_mat) @ refwave + + tcircuit, phase = trotterize(qubit_hamiltonian, trotter_order=1, n_trotter_steps=1, time=time, + return_phase=True) + _, wavefunc = sim.simulate(tcircuit, return_statevector=True, initial_statevector=refwave) + wavefunc *= phase + overlap = np.dot(np.conj(evolve_exact), wavefunc) + self.assertAlmostEqual(overlap, 1.0, delta=1e-3) + + def test_trotterize_different_order_and_steps(self): + """ Verify that the time evolution is correct for different orders and number of steps + with a qubit_hamiltonian input""" + + time = 0.2 + mapping = "bk" + reference_circuit = get_reference_circuit(n_spinorbitals=mol_H4_sto3g.n_active_sos, + n_electrons=mol_H4_sto3g.n_active_electrons, + mapping=mapping, + up_then_down=True) + _, refwave = sim.simulate(reference_circuit, return_statevector=True) + + qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=fermion_operator, + mapping=mapping, + n_spinorbitals=mol_H4_sto3g.n_active_sos, + n_electrons=mol_H4_sto3g.n_active_electrons, + up_then_down=True) + + ham_mat = get_sparse_operator(qubit_hamiltonian).toarray() + evolve_exact = expm(-1j * time * ham_mat) @ refwave + + for trotter_order, n_trotter_steps in [(1, 1), (2, 1), (1, 2)]: + + tcircuit, phase = trotterize(qubit_hamiltonian, time, n_trotter_steps, trotter_order, return_phase=True) + _, wavefunc = sim.simulate(tcircuit, return_statevector=True, initial_statevector=refwave) + wavefunc *= phase + overlap = np.dot(np.conj(evolve_exact), wavefunc) + self.assertAlmostEqual(overlap, 1.0, delta=1e-3) + + def test_trotterize_fermionic_input_different_times(self): + """ Verify that the time evolution is correct for a FermionOperator input with different times + for each term + """ + + mapping = "jw" + # generate Hermitian FermionOperator + fermion_operators = [FermionOperator("0^ 3", 0.5) + FermionOperator("3^ 0", 0.5), + FermionOperator("1^ 2", 0.5) + FermionOperator("2^ 1", 0.5), + FermionOperator("1^ 3", 0.5) + FermionOperator("3^ 1", 0.5)] + + # time is twice as long as each Hermitian Operator has two terms + time = {((0, 1), (3, 0)): 0.1, ((3, 1), (0, 0)): 0.1, ((1, 1), (2, 0)): 0.2, ((2, 1), (1, 0)): 0.2, + ((1, 1), (3, 0)): 0.3, ((3, 1), (1, 0)): 0.3} + + # Build referenc circuit and obtain reference wavefunction + reference_circuit = Circuit([Gate("X", 0), Gate("X", 3)]) + _, refwave = sim.simulate(reference_circuit, return_statevector=True) + + evolve_exact = refwave + total_fermion_operator = FermionOperator() + # evolve each term separately and apply to resulting wavefunction + for i in range(3): + total_fermion_operator += fermion_operators[i] + qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=fermion_operators[i], + mapping=mapping) + ham_mat = get_sparse_operator(qubit_hamiltonian, n_qubits=4).toarray() + evolve_exact = expm(-1j * time[next(iter(fermion_operators[i].terms))] * ham_mat) @ evolve_exact + + # Apply trotter-suzuki steps using different times for each term + tcircuit, phase = trotterize(total_fermion_operator, trotter_order=1, n_trotter_steps=1, time=time, return_phase=True) + _, wavefunc = sim.simulate(tcircuit, return_statevector=True, initial_statevector=refwave) + wavefunc *= phase + + overlap = np.dot(np.conj(evolve_exact), wavefunc) + self.assertAlmostEqual(overlap, 1.0, delta=1e-3) + + def test_trotterize_qubit_input_different_times(self): + """ Verify that the time evolution is correct for a QubitOperator input with different times + for each term + """ + + qubit_operator_list = [QubitOperator("X0 Y1", 0.5), QubitOperator("Y1 Z2", 0.5), QubitOperator("Y2 X3", 0.5)] + + time = {((0, 'X'), (1, 'Y')): 0.1, ((1, 'Y'), (2, 'Z')): 0.2, ((2, 'Y'), (3, 'X')): 0.3} + + # Generate initial wavefunction + reference_circuit = Circuit([Gate("X", 0), Gate("X", 3)]) + _, refwave = sim.simulate(reference_circuit, return_statevector=True) + + # Exactly evolve for each time step + evolve_exact = refwave + for i in range(3): + ham_mat = get_sparse_operator(qubit_operator_list[i], n_qubits=4).toarray() + evolve_exact = expm(-1j * time[next(iter(qubit_operator_list[i].terms))] * ham_mat) @ evolve_exact + + # Apply trotter-suzuki with different times for each qubit operator term + total_qubit_operator = QubitOperator() + for qu_op in reversed(qubit_operator_list): + total_qubit_operator += qu_op + tcircuit, phase = trotterize(total_qubit_operator, trotter_order=2, n_trotter_steps=2, time=time, return_phase=True) + _, wavefunc = sim.simulate(tcircuit, return_statevector=True, initial_statevector=refwave) + wavefunc *= phase + overlap = np.dot(np.conj(evolve_exact), wavefunc) + self.assertAlmostEqual(overlap, 1.0, delta=1e-3) + + def test_qft_by_phase_estimation(self): + """Test get_qft_circuit by applying phase-estimation to a 1-qubit operator with eigenvalue -1, i.e. phi=1/2""" + n_qubits = 4 + qubit_list = [2, 1, 0] + # Generate state with eigenvalue -1 of X operator exp(2*pi*i*phi) phi=1/2 + gate_list = [Gate("X", target=n_qubits-1), Gate("H", target=n_qubits-1)] + + # Generate phase-estimation circuit with three registers + pe_circuit = Circuit(gate_list, n_qubits=n_qubits) + qft = get_qft_circuit(qubit_list, n_qubits=n_qubits) + pe_circuit += qft + controlled_unitaries = [] + for i, qubit in enumerate(qubit_list): + for j in range(2**i): + controlled_unitaries += [Gate("CNOT", target=n_qubits-1, control=qubit)] + pe_circuit += Circuit(controlled_unitaries, n_qubits=n_qubits) + iqft = get_qft_circuit(qubit_list, n_qubits=n_qubits, inverse=True) + pe_circuit += iqft + + # simulate starting state frequency is {"0001": 0.5, "0000": 0.5} + freqs, _ = sim.simulate(pe_circuit) + # phase is added to first three qubits with value 100 = 1 * 1/2 + 0 * 1/4 + 0 * 1/8 + # while keeping last qubit unchanged. Therefore, the target frequency dictionary is target_freq_dict + target_freq_dict = {"1000": 0.5, "1001": 0.5} + assert_freq_dict_almost_equal(target_freq_dict, freqs, atol=1.e-7) + + def test_controlled_time_evolution_by_phase_estimation(self): + """ Verify that the time evolution is correct for a QubitOperator input with different times + for each term + """ + + # Generate qubit operator with state 9 having eigenvalue 0.25 + qu_op = (QubitOperator("X0 X1", 0.125) + QubitOperator("Y1 Y2", 0.125) + QubitOperator("Z2 Z3", 0.125) + + QubitOperator("", 0.125)) + + ham_mat = get_sparse_operator(qu_op).toarray() + _, wavefunction = eigh(ham_mat) + + # Append four qubits in the zero state to eigenvector 9 + wave_9 = wavefunction[:, 9] + for i in range(4): + wave_9 = np.kron(wave_9, np.array([1, 0])) + + n_qubits = 8 + + qubit_list = [7, 6, 5, 4] + + qft = get_qft_circuit(qubit_list, n_qubits=n_qubits) + pe_circuit = qft + for i, qubit in enumerate(qubit_list): + u_circuit = trotterize(qu_op, trotter_order=1, n_trotter_steps=10, time=-2*np.pi, control=qubit) + for j in range(2**i): + pe_circuit += u_circuit + iqft = get_qft_circuit(qubit_list, n_qubits=n_qubits, inverse=True) + pe_circuit += iqft + + freqs, _ = sim.simulate(pe_circuit, initial_statevector=wave_9) + + # Trace out first 4 dictionary amplitudes, only care about final 4 indices + trace_freq = dict() + for key, value in freqs.items(): + trace_freq[key[-4:]] = trace_freq.get(key[-4:], 0) + value + + # State 9 has eigenvalue 0.25 so return should be 0100 (0*1/2 + 1*1/4 + 0*1/8 + 0*1/16) + self.assertAlmostEqual(trace_freq["0100"], 1.0, delta=2) + + def test_controlled_swap(self): + cswap_circuits = [Circuit([Gate("CSWAP", target=[1, 2], control=0)]), + Circuit(controlled_swap_to_XX_gates(0, 1, 2))] + + for cswap_circuit in cswap_circuits: + # initialize in "110", returns "101" + init_gates = [Gate("X", target=0), Gate("X", target=1)] + circuit = Circuit(init_gates, n_qubits=3) + cswap_circuit + freqs, _ = sim.simulate(circuit, return_statevector=True) + assert_freq_dict_almost_equal({"101": 1.0}, freqs, atol=1.e-7) + + # initialize in "010" returns "010" + init_gates = [Gate("X", target=1)] + circuit = Circuit(init_gates, n_qubits=3) + cswap_circuit + freqs, _ = sim.simulate(circuit, return_statevector=True) + assert_freq_dict_almost_equal({"010": 1.0}, freqs, atol=1.e-7) + + def test_derangement_circuit_by_estimating_pauli_string(self): + """ Verify that tr(rho^3 pa) for a pauliword pa is correct. + Uses the exponential error suppression circuit + """ + + qu_op = QubitOperator("X0 Y1", 0.125) + QubitOperator("Y1 Y2", 0.125) + QubitOperator("Z2 Z3", 0.125) + pa = QubitOperator("X0 X1 X2 X3", 1) + + ham_mat = get_sparse_operator(qu_op).toarray() + _, wavefunction = eigh(ham_mat) + pamat = get_sparse_operator(pa).toarray() + + mixed_wave = np.sqrt(3)/2*wavefunction[:, -1] + 1/2*wavefunction[:, 0] + mixed_wave_3 = np.kron(np.kron(mixed_wave, mixed_wave), mixed_wave) + full_start_vec = np.kron(mixed_wave_3, np.array([1, 0])) + rho = np.outer(mixed_wave, mixed_wave) + rho3 = rho @ rho @ rho + exact = np.trace(rho3 @ pamat) + + n_qubits = 13 + + qubit_list = [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]] + rho3_pa_circuit = Circuit([Gate("H", target=12)], n_qubits=n_qubits) + derange_circuit = derangement_circuit(qubit_list, control=12, n_qubits=n_qubits) + cpas = controlled_pauliwords(pa, control=12, n_qubits=n_qubits) + rho3_pa_circuit += derange_circuit + cpas[0] + rho3_pa_circuit += Circuit([Gate("H", target=12)], n_qubits=n_qubits) + + exp_op = QubitOperator("Z12", 1) + measured = sim.get_expectation_value(exp_op, rho3_pa_circuit, initial_statevector=full_start_vec) + self.assertAlmostEqual(measured, exact, places=6) + + +if __name__ == "__main__": + unittest.main() diff --git a/tangelo/toolboxes/ansatz_generator/uccsd.py b/tangelo/toolboxes/ansatz_generator/uccsd.py index 8e5c26d8c..58c604943 100644 --- a/tangelo/toolboxes/ansatz_generator/uccsd.py +++ b/tangelo/toolboxes/ansatz_generator/uccsd.py @@ -36,7 +36,7 @@ from tangelo.linq import Circuit from .ansatz import Ansatz -from .ansatz_utils import pauliword_to_circuit +from .ansatz_utils import exp_pauliword_to_gates from ._unitary_cc_openshell import uccsd_openshell_paramsize, uccsd_openshell_generator from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit @@ -166,7 +166,7 @@ def build_circuit(self, var_params=None): pauli_words_gates = [] self.pauli_to_angles_mapping = dict() for i, (pauli_word, coef) in enumerate(pauli_words): - pauli_words_gates += pauliword_to_circuit(pauli_word, coef) + pauli_words_gates += exp_pauliword_to_gates(pauli_word, coef) self.pauli_to_angles_mapping[pauli_word] = i uccsd_circuit = Circuit(pauli_words_gates) diff --git a/tangelo/toolboxes/ansatz_generator/upccgsd.py b/tangelo/toolboxes/ansatz_generator/upccgsd.py index b78093777..c4c2879b6 100644 --- a/tangelo/toolboxes/ansatz_generator/upccgsd.py +++ b/tangelo/toolboxes/ansatz_generator/upccgsd.py @@ -29,7 +29,7 @@ from tangelo.linq import Circuit from .ansatz import Ansatz -from .ansatz_utils import pauliword_to_circuit +from .ansatz_utils import exp_pauliword_to_gates from ._unitary_cc_paired import get_upccgsd from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit @@ -156,7 +156,7 @@ def build_circuit(self, var_params=None): # Obtain quantum circuit through trivial trotterization of the qubit operator for each current_k for i, (pauli_word, coef) in enumerate(pauli_words): - pauli_words_gates += pauliword_to_circuit(pauli_word, coef) + pauli_words_gates += exp_pauliword_to_gates(pauli_word, coef) self.pauli_to_angles_mapping[current_k][pauli_word] = i + sum_prev_qubit_terms[current_k] sum_prev_qubit_terms[current_k + 1] = len(qubit_op.terms.items())