diff --git a/CHANGELOG.md b/CHANGELOG.md index 192b147dc..2c94fb5e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,8 +12,6 @@ This file documents the main changes between versions of the code. - iQCC ansatz for VQE - IonQConnection class and notebook, to facilitate experiments through IonQ's API - FCISolver active space selection / frozen orbitals: restrictions for half-empty orbitals -- HybridOperator for speedup for QubitOperator on certain operations in stabilizer notation -- Support for symmetry in pyscf computations ### Changed diff --git a/tangelo/algorithms/variational/__init__.py b/tangelo/algorithms/variational/__init__.py index 9d46fddcc..e14ef27b4 100644 --- a/tangelo/algorithms/variational/__init__.py +++ b/tangelo/algorithms/variational/__init__.py @@ -17,3 +17,4 @@ from .sa_vqe_solver import SA_VQESolver from .sa_oo_vqe_solver import SA_OO_Solver from .iqcc_solver import iQCC_solver +from .iqcc_ilc_solver import iQCC_ILC_solver diff --git a/tangelo/algorithms/variational/iqcc_ilc_solver.py b/tangelo/algorithms/variational/iqcc_ilc_solver.py new file mode 100644 index 000000000..952f7c1b3 --- /dev/null +++ b/tangelo/algorithms/variational/iqcc_ilc_solver.py @@ -0,0 +1,302 @@ +# 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 implements the iQCC-ILC VQE procedure of Refs. 1-2. It is +an iterative and variational approach that combines an ansatz defined +as an exponentiated involutory linear combination (ILC) of mutually +anticommuting Pauli word generators with the QCC ansatz. A small number +of iterations are performed with the ILC ansatz prior to a single energy +evaluation with the QCC ansatz. The advantage of this method over the +iQCC VQE procedure is that Hamiltonian dressing after each iteration +with the set of ILC generators results in quadratic growth of the +number of terms, which is an improvement over the exponential growth +that occurs when QCC generators are used. +Refs: + 1. R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov. + arXiv:2002.05701v1, 2020, 1–10. + 2. R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov. + J. Chem. Theory Comput. 2021, 17, 1, 66–78. +""" + +from tangelo.linq import Simulator +from tangelo.toolboxes.ansatz_generator.ilc import ILC +from tangelo.toolboxes.ansatz_generator.qcc import QCC +from tangelo.algorithms.variational.vqe_solver import VQESolver +from tangelo.toolboxes.ansatz_generator._qubit_ilc import ilc_op_dress + + +class iQCC_ILC_solver: + """The iQCC-ILC-VQE solver class combines the both the ILC and ILC ansatze + Classes with the VQESolver class to perform an iterative and variational + procedure to compute the total iQCC-ILC energy for a given Hamiltonian. + The algorithm is outlined below: + (1) For a user-specified number of iterations, compute the ILC energy: + (a) prepare/purify the QMF wave function, obtain the ACS of ILC + generators, and initialize the ILC parameter set; + (b) simulate the ILC energy through VQE minimization + (c) dress the qubit Hamiltonian with the set of ILC generators and + optimal parameters; optional: compress the dressed Hamiltonian + via a technique using the Frobenius norm + (2) With the ILC dressed Hamiltonian, obtain the DIS of QCC generators, + and initialize QCC parameters + (3) Perform a single VQE minimization of the QCC energy functional to + obtain the final iQCC-ILC energy. + Attributes: + molecule (SecondQuantizedMolecule): The molecular system. + qubit_mapping (str): One of the supported qubit mapping identifiers. Default, "jw". + up_then_down (bool): Change basis ordering putting all spin up orbitals first, + followed by all spin down. Default, False. + initial_var_params (str or array-like): Initial values of the variational parameters + for the classical optimizer. + backend_options (dict): Parameters to build the tangelo.linq Simulator + class. + penalty_terms (dict): Parameters for penalty terms to append to target + qubit Hamiltonian (see penaly_terms for more details). + ilc_ansatz_options (dict): Parameters for ILC ansatz (see ILC ansatz + file for details). + qcc_ansatz_options (dict): Parameters for QCC ansatz (see QCC ansatz + file for details). + qubit_hamiltonian (QubitOperator-like): Self-explanatory. + max_ilc_iter (int): maximum number of ILC iterations allowed before termination. + Default, 3. + compress_qubit_ham (bool): controls whether the qubit Hamiltonian is compressed + after dressing with the current set of generators at the end of each ILC iteration. + Default, False. + compress_eps (float): parameter required for compressing intermediate ILC Hamiltonians + using the Froebenius norm. Discarding terms in this manner will not alter the + eigenspectrum of intermediate Hamiltonians by more than compress_eps. + Default, 1.59e-3 Hartree. + verbose (bool): Flag for verbosity. Default, False. + """ + + def __init__(self, opt_dict): + + default_backend_options = {"target": None, "n_shots": None, "noise_model": None} + default_options = {"molecule": None, + "qubit_mapping": "jw", + "up_then_down": False, + "initial_var_params": None, + "backend_options": default_backend_options, + "penalty_terms": None, + "ilc_ansatz_options": dict(), + "qcc_ansatz_options": dict(), + "qubit_hamiltonian": None, + "max_ilc_iter": 3, + "compress_qubit_ham": False, + "compress_eps": 1.59e-3, + "verbose": False} + + # Initialize with default values + self.__dict__ = default_options + # Overwrite default values with user-provided ones, if they correspond to a valid keyword + for param, val in opt_dict.items(): + if param in default_options: + setattr(self, param, val) + else: + raise KeyError(f"The keyword {param} is not available in self.__class__.__name__.") + + if not self.molecule and not self.qubit_hamiltonian: + raise ValueError("An instance of SecondQuantizedMolecule or QubitOperator " + "is required for initializing self.__class__.__name__.") + + # initialize variables and lists to store useful data from each ILC-VQE iteration + self.energies = [] + self.iteration = 0 + self.terminate_ilc = False + self.qmf_energy = None + self.ilc_ansatz = None + self.qcc_ansatz = None + self.vqe_solver = None + self.vqe_solver_options = None + self.final_optimal_energy = None + self.final_optimal_qmf_params = None + self.final_optimal_ilc_params = None + self.final_optimal_qcc_params = None + + def build(self): + """Builds the underlying objects required to run the ILC-VQE algorithm.""" + + # instantiate the ILC ansatz but do not build it here because vqe_solver builds it + self.ilc_ansatz = ILC(self.molecule, self.qubit_mapping, self.up_then_down, **self.ilc_ansatz_options) + + # build an instance of VQESolver with options that remain fixed during the ILC-VQE routine + self.vqe_solver_options = {"qubit_hamiltonian": self.ilc_ansatz.qubit_ham, + "qubit_mapping": self.qubit_mapping, + "ansatz": self.ilc_ansatz, + "initial_var_params": self.initial_var_params, + "backend_options": self.backend_options, + "penalty_terms": self.penalty_terms, + "up_then_down": self.up_then_down, + "verbose": self.verbose} + + self.vqe_solver = VQESolver(self.vqe_solver_options) + self.vqe_solver.build() + + def simulate(self): + """Executes the ILC-VQE algorithm. During each iteration, a ILC-VQE minimization + is performed with the current set of generators, amplitudes, and qubit Hamiltonian.""" + + # initialize quantities and compute the QMF energy + sim = Simulator() + self.qmf_energy = sim.get_expectation_value(self.ilc_ansatz.qubit_ham, self.ilc_ansatz.qmf_circuit) + if self.verbose: + print(f"The qubit mean field energy = {self.qmf_energy}") + + # perform self.max_ilc_iter ILC-VQE minimizations; + e_ilc = 0. + while not self.terminate_ilc: + # check that the ACS has at least one ILC generator to use; if not, terminate + if self.ilc_ansatz.acs and self.ilc_ansatz.var_params.any(): + e_ilc = self.vqe_solver.simulate() + else: + self.terminate_ilc = True + if self.verbose: + print("Terminating the ILC-VQE solver: the ACS of ILC generators is empty.") + # update ILC-VQE simulation data + if not self.terminate_ilc: + self._update_ilc_solver(e_ilc) + + # perform a single QCC-VQE minimization to obtain the final iQCC-ILC energy + # need to rebuild VQE Solver for the QCC ansatz first + self._build_qcc() + + # check that the DIS has at least one QCC generator to use + if self.qcc_ansatz.dis and self.qcc_ansatz.var_params.any(): + e_iqcc_ilc = self.vqe_solver.simulate() + self._update_qcc_solver(e_iqcc_ilc) + else: + if self.verbose: + print("Terminating the iQCC-ILC solver without evaluating the " + "the final QCC energy because the DIS of QCC generators " + "is empty for the given Hamiltonian.") + + return self.energies[-1] + + def get_resources(self): + """Returns the quantum resource estimates for the final + ILC-QCC-VQE iteration.""" + + return self.vqe_solver.get_resources() + + def _build_qcc(self): + """Builds the underlying QCC objects required to run the second part of the + iQCC-ILC-VQE algorithm.""" + + # instantiate the QCC ansatz with the final ILC dressed Hamiltonian and optimal QMF params + self.qcc_ansatz_options["qmf_var_params"] = self.final_optimal_qmf_params + self.qcc_ansatz_options["qubit_ham"] = self.ilc_ansatz.qubit_ham + self.qcc_ansatz = QCC(self.molecule, self.qubit_mapping, self.up_then_down, **self.qcc_ansatz_options) + + # build an instance of VQESolver with options that remain fixed during the ILC-VQE routine + self.vqe_solver_options = {"qubit_mapping": self.qubit_mapping, + "ansatz": self.qcc_ansatz, + "initial_var_params": self.initial_var_params, + "backend_options": self.backend_options, + "penalty_terms": self.penalty_terms, + "up_then_down": self.up_then_down, + "qubit_hamiltonian": self.ilc_ansatz.qubit_ham, + "verbose": self.verbose} + + self.vqe_solver = VQESolver(self.vqe_solver_options) + self.vqe_solver.build() + + def _update_ilc_solver(self, e_ilc): + """This function serves several purposes for the ILC-VQE solver + part of the iQCC-ILC algorithm: + + (1) Updates the ILC energy, generators, QMF Bloch angles, + ILC amplitudes, circuits, number of qubit Hamiltonian terms, + and quantum resource estimates; + (2) Dresses and compresses (optional) the qubit Hamiltonian + with the generators and optimal amplitudes for the current + iteration; + (3) Prepares for the next iteration by rebuilding the DIS & ACS, + reinitializing ILC parameters for the new set of generators, + generating the circuit, and updating the classical optimizer. + + Args: + e_ilc (float): the total ILC ansatz energy at the current iteration. + """ + + # get the optimal variational parameters and split them for qmf and ilc + n_qubits = self.ilc_ansatz.n_qubits + optimal_qmf_var_params = self.vqe_solver.optimal_var_params[:2*n_qubits] + optimal_ilc_var_params = self.vqe_solver.optimal_var_params[2*n_qubits:] + + # update energy list and iteration number + self.energies.append(e_ilc) + self.iteration += 1 + + if self.verbose: + print(f"Iteration # {self.iteration}") + print(f"ILC total energy = {e_ilc} Eh") + print(f"ILC correlation energy = {e_ilc-self.qmf_energy} Eh") + print(f"Optimal QMF variational parameters = {optimal_qmf_var_params}") + print(f"Optimal ILC variational parameters = {optimal_ilc_var_params}") + print(f"# of ILC generators = {len(self.ilc_ansatz.acs)}") + print(f"ILC generators = {self.ilc_ansatz.acs}") + print(f"ILC resource estimates = {self.get_resources()}") + if self.iteration == 1: + n_qham_terms = self.ilc_ansatz.qubit_ham.get_max_number_hamiltonian_terms(n_qubits) + print(f"Maximum number of qubit Hamiltonian terms = {n_qham_terms}") + + # dress and (optionally) compress the qubit Hamiltonian for the next iteration + self.ilc_ansatz.qubit_ham = ilc_op_dress(self.ilc_ansatz.qubit_ham, self.ilc_ansatz.acs, + optimal_ilc_var_params) + if self.compress_qubit_ham: + self.ilc_ansatz.qubit_ham.frobenius_norm_compression(self.compress_eps, n_qubits) + + # set dis, acs, and var_params to none to rebuild the dis & acs and initialize new parameters + self.ilc_ansatz.dis = None + self.ilc_ansatz.acs = None + self.ilc_ansatz.var_params = None + self.ilc_ansatz.build_circuit() + + self.vqe_solver.qubit_hamiltonian = self.ilc_ansatz.qubit_ham + self.vqe_solver.initial_var_params = self.ilc_ansatz.var_params + + if self.iteration == self.max_ilc_iter: + self.terminate_ilc = True + self.final_optimal_qmf_params = optimal_qmf_var_params + self.final_optimal_ilc_params = optimal_ilc_var_params + if self.verbose: + print(f"Terminating the ILC-VQE solver after {self.max_ilc_iter} iterations.") + + def _update_qcc_solver(self, e_iqcc_ilc): + """Updates after the final QCC energy evaluation with the final ILC dressed + Hamiltonian. + + Args: + e_iqcc_ilc (float): the final iQCC-ILC ansatz energy. + """ + + # get the optimal variational parameters and split them for qmf and qcc ansatze + n_qubits = self.qcc_ansatz.n_qubits + self.final_optimal_qmf_var_params = self.vqe_solver.optimal_var_params[:2*n_qubits] + self.final_optimal_qcc_var_params = self.vqe_solver.optimal_var_params[2*n_qubits:] + + # update energy list + self.final_optimal_energy = e_iqcc_ilc + self.energies.append(e_iqcc_ilc) + + if self.verbose: + print("Final iQCC-ILC VQE simulation.") + print(f"iQCC-ILC total energy = {e_iqcc_ilc} Eh") + print(f"iQCC-ILC correlation energy = {e_iqcc_ilc-self.qmf_energy} Eh") + print(f"Optimal QMF variational parameters = {self.final_optimal_qmf_var_params}") + print(f"Optimal QCC variational parameters = {self.final_optimal_qcc_var_params}") + print(f"# of QCC generators = {len(self.qcc_ansatz.dis)}") + print(f"QCC generators = {self.qcc_ansatz.dis}") + print(f"QCC resource estimates = {self.get_resources()}") diff --git a/tangelo/algorithms/variational/iqcc_solver.py b/tangelo/algorithms/variational/iqcc_solver.py index cd90459e3..ac30cecd2 100644 --- a/tangelo/algorithms/variational/iqcc_solver.py +++ b/tangelo/algorithms/variational/iqcc_solver.py @@ -21,9 +21,7 @@ for the iQCC approach relative to the native QCC method. A caveat is that after each iteration, the qubit Hamiltonian is dressed with the generators and optimal parameters, the result of which is an -exponential growth of the number of terms. A technique also described -in Ref. 1 can be utilized to address this issue by discarding some -terms based on the Frobenius norm of the Hamiltonian. +exponential growth of the number of terms. Refs: 1. I. G. Ryabinkin, R. A. Lang, S. N. Genin, and A. F. Izmaylov. @@ -88,7 +86,7 @@ class iQCC_solver: using the Froebenius norm. Discarding terms in this manner will not alter the eigenspeectrum of intermediate Hamiltonians by more than compress_eps. Default, 1.59e-3 Hartree. - verbose (bool): Flag for verbosity of iQCCsolver. Default, False. + verbose (bool): Flag for verbosity. Default, False. """ def __init__(self, opt_dict): @@ -116,10 +114,10 @@ def __init__(self, opt_dict): if param in default_options: setattr(self, param, val) else: - raise KeyError(f"Keyword :: {param}, not available in iQCCsolver") + raise KeyError(f"Keyword {param} not available in self.__class__.__name__.") if not self.molecule: - raise ValueError("An instance of SecondQuantizedMolecule is required for initializing iQCCsolver.") + raise ValueError("An instance of SecondQuantizedMolecule is required for initializing self.__class__.__name__.") # initialize variables and lists to store useful data from each iQCC-VQE iteration self.energies = [] @@ -149,6 +147,7 @@ def build(self): "up_then_down": self.up_then_down, "qubit_hamiltonian": self.qubit_hamiltonian, "verbose": self.verbose} + self.vqe_solver = VQESolver(self.vqe_solver_options) self.vqe_solver.build() @@ -233,8 +232,19 @@ def _update_iqcc_solver(self, delta_eqcc): optimal_qmf_var_params = self.vqe_solver.optimal_var_params[:2*n_qubits] optimal_qcc_var_params = self.vqe_solver.optimal_var_params[2*n_qubits:] - # update all lists with data from the current iteration + # update energy list and iteration number self.energies.append(self.vqe_solver.optimal_energy) + self.iteration += 1 + + if self.verbose: + print(f"Iteration # {self.iteration}") + print(f"iQCC total energy = {self.vqe_solver.optimal_energy} Eh") + print(f"iQCC correlation energy = {self.vqe_solver.optimal_energy-self.qmf_energy} Eh") + print(f"Optimal QMF variational parameters = {optimal_qmf_var_params}") + print(f"Optimal QCC variational parameters = {optimal_qcc_var_params}") + print(f"Number of iQCC generators = {len(self.qcc_ansatz.dis)}") + print(f"iQCC generators = {self.qcc_ansatz.dis}") + print(f"iQCC resource estimates = {self.get_resources()}") # dress and (optionally) compress the qubit Hamiltonian self.qcc_ansatz.qubit_ham = qcc_op_dress(self.qcc_ansatz.qubit_ham, self.qcc_ansatz.dis, @@ -248,18 +258,6 @@ def _update_iqcc_solver(self, delta_eqcc): self.qcc_ansatz.build_circuit() self.vqe_solver.initial_var_params = self.qcc_ansatz.var_params - self.iteration += 1 - - if self.verbose: - print(f"Iteration # {self.iteration}") - print(f"iQCC total energy = {self.vqe_solver.optimal_energy} Eh") - print(f"iQCC correlation energy = {self.vqe_solver.optimal_energy-self.qmf_energy} Eh") - print(f"Optimal QMF variational parameters = {optimal_qmf_var_params}") - print(f"Optimal QCC variational parameters = {optimal_qcc_var_params}") - print(f"Number of iQCC generators = {len(self.qcc_ansatz.dis)}") - print(f"iQCC generators = {self.qcc_ansatz.dis}") - print(f"iQCC resource estimates = {self.get_resources()}") - if abs(delta_eqcc) < self.deqcc_thresh or self.iteration == self.max_iqcc_iter: self.converged = True self.final_optimal_energy = self.vqe_solver.optimal_energy diff --git a/tangelo/algorithms/variational/tests/test_iqcc_ilc_solver.py b/tangelo/algorithms/variational/tests/test_iqcc_ilc_solver.py new file mode 100644 index 000000000..6c3fad596 --- /dev/null +++ b/tangelo/algorithms/variational/tests/test_iqcc_ilc_solver.py @@ -0,0 +1,91 @@ +# 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. + +"""Unit tests for the closed-shell and restricted open-shell iQCC-ILC solver.""" + +import unittest + +from tangelo.algorithms.variational import iQCC_ILC_solver +from tangelo.molecule_library import mol_H2_sto3g, mol_H4_sto3g_symm, mol_H4_cation_sto3g + + +class iQCC_ILC_solver_test(unittest.TestCase): + """Unit tests for the iQCC_ILC_solver class. Examples for both closed-shell + and restricted open-shell iQCC-ILC are provided via H4 and H4+. + """ + + @staticmethod + def test_build_success(): + """Test instantation of the iQCC-ILC solver with user-defined input.""" + + iqcc_ilc_options = {"molecule": mol_H2_sto3g, + "qubit_mapping": "scbk", + "up_then_down": True, + "max_ilc_iter": 25, + "compress_qubit_ham": True, + "compress_eps": 1e-4} + + iqcc_ilc = iQCC_ILC_solver(iqcc_ilc_options) + iqcc_ilc.build() + + def test_build_fail(self): + """Test that instantation of the iQCC-ILC solver fails without input of a molecule.""" + + iqcc_ilc_options = {"max_ilc_iter": 15} + self.assertRaises(ValueError, iQCC_ILC_solver, iqcc_ilc_options) + + def test_iqcc_ilc_h4(self): + """Test the energy after 1 iteration for H4.""" + + ilc_ansatz_options = {"max_ilc_gens": None} + qcc_ansatz_options = {"max_qcc_gens": None} + + iqcc_ilc_options = {"molecule": mol_H4_sto3g_symm, + "qubit_mapping": "scbk", + "up_then_down": True, + "ilc_ansatz_options": ilc_ansatz_options, + "qcc_ansatz_options": qcc_ansatz_options, + "max_ilc_iter": 1} + + iqcc_ilc_solver = iQCC_ILC_solver(iqcc_ilc_options) + iqcc_ilc_solver.build() + iqcc_ilc_energy = iqcc_ilc_solver.simulate() + + self.assertAlmostEqual(iqcc_ilc_energy, -1.976817, places=4) + + def test_iqcc_ilc_h4_cation(self): + """Test the energy after 2 iterations for H4+ using the maximum + number of generators and compressing the qubit Hamiltonian.""" + + ilc_ansatz_options = {"max_ilc_gens": None} + qcc_ansatz_options = {"max_qcc_gens": None} + + iqcc_ilc_options = {"molecule": mol_H4_cation_sto3g, + "qubit_mapping": "scbk", + "up_then_down": True, + "ilc_ansatz_options": ilc_ansatz_options, + "qcc_ansatz_options": qcc_ansatz_options, + "max_ilc_iter": 2, + "compress_qubit_ham": True, + "compress_eps": 1e-4} + + iqcc_ilc_solver = iQCC_ILC_solver(iqcc_ilc_options) + iqcc_ilc_solver.build() + iqcc_ilc_energy = iqcc_ilc_solver.simulate() + + self.assertAlmostEqual(iqcc_ilc_energy, -1.638197, places=4) + + +if __name__ == "__main__": + unittest.main() diff --git a/tangelo/algorithms/variational/tests/test_iqcc_solver.py b/tangelo/algorithms/variational/tests/test_iqcc_solver.py index 2c18693ef..cb1242786 100644 --- a/tangelo/algorithms/variational/tests/test_iqcc_solver.py +++ b/tangelo/algorithms/variational/tests/test_iqcc_solver.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -"""Unit tests for the closed-shell and restricted open-shell iQCC-VQE Solver. """ +"""Unit tests for the closed-shell and restricted open-shell iQCC VQE Solver. """ import unittest diff --git a/tangelo/toolboxes/ansatz_generator/_qubit_cc.py b/tangelo/toolboxes/ansatz_generator/_qubit_cc.py index ee6f04f27..4a0d5d74f 100644 --- a/tangelo/toolboxes/ansatz_generator/_qubit_cc.py +++ b/tangelo/toolboxes/ansatz_generator/_qubit_cc.py @@ -161,17 +161,16 @@ def get_gens_from_idxs(group_idxs): return dis_group_gens -def build_qcc_qubit_op(dis_gens, amplitudes): - """Returns the QCC operator by selecting n_var_params generators from the DIS. +def build_qcc_qubit_op(dis_gens, taus): + """Returns the QCC operator with generators selected from the DIS. The QCC operator is constructed as a linear combination of generators using the - parameter set {tau} as coefficients: QCC operator = -0.5 * SUM_k P_k * tau_k. - The exponentiated QCC operator, U = PROD_k exp(-0.5j * tau_k * P_k), is used to - build the circuit. + taus parameter set: QCC operator = -0.5 * SUM_k P_k * tau_k. The exponentiated + QCC operator, U = PROD_k exp(-0.5j * tau_k * P_k), is used to build the circuit. Args: dis_gens (list of QubitOperator): The list of QCC Pauli word generators selected from a user-specified number of characteristic DIS groups. - amplitudes (list or numpy array of float): The QCC variational parameters + taus (list or numpy array of float): The QCC variational parameters arranged such that their ordering matches the order of dis_gens. Returns: @@ -180,25 +179,23 @@ def build_qcc_qubit_op(dis_gens, amplitudes): qubit_op = QubitOperator.zero() for i, dis_gen in enumerate(dis_gens): - qubit_op -= 0.5 * amplitudes[i] * dis_gen + qubit_op -= 0.5 * taus[i] * dis_gen qubit_op.compress() return qubit_op -def qcc_op_dress(qubit_op, dis_gens, amplitudes): +def qcc_op_dress(qubit_op, dis_gens, taus): """Performs canonical transformation of a qubit operator with the set of QCC - generators and amplitudes for the current iteration. For an operator with M terms - each transformation results in exponential growth of the number terms. This growth - can be approximated as M * (3 / 2) ^ n_g, where n_g is the number of QCC generators + generators and amplitudes. For an operator with M terms, each transformation + results in exponential growth of the number terms. This growth can be + approximated as M * (3 / 2) ^ n_g, where n_g is the number of QCC generators selected for the ansatz at the current iteration. Args: - qubit_op (QubitOperator): A qubit operator (e.g., a molecular Hamiltonian or the - electronic spin and number operators) that was previously dressed by canonical - transformation with the QCC generators and amplitudes at the current iteration. + qubit_op (QubitOperator): A qubit operator to be dressed. dis_gens (list of QubitOperator): The list of QCC Pauli word generators - selected from a user-specified number of characteristic DIS groups. - amplitudes (list or numpy array of float): The QCC variational parameters + selected from characteristic DIS groups. + taus (list or numpy array of float): The QCC variational parameters arranged such that their ordering matches the ordering of dis_gens. Returns: @@ -206,7 +203,6 @@ def qcc_op_dress(qubit_op, dis_gens, amplitudes): """ for i, gen in enumerate(dis_gens): - comm = commutator(qubit_op, gen) - qubit_op += .5 * ((1. - cos(amplitudes[i])) * gen - 1j * sin(amplitudes[i])) * comm + qubit_op += .5 * ((1. - cos(taus[i])) * gen - 1j * sin(taus[i])) * commutator(qubit_op, gen) qubit_op.compress() return qubit_op diff --git a/tangelo/toolboxes/ansatz_generator/_qubit_ilc.py b/tangelo/toolboxes/ansatz_generator/_qubit_ilc.py index 419c17cc2..f5d5fed92 100644 --- a/tangelo/toolboxes/ansatz_generator/_qubit_ilc.py +++ b/tangelo/toolboxes/ansatz_generator/_qubit_ilc.py @@ -22,35 +22,37 @@ arXiv:2002.05701v1, 2020, 1–10. 2. R. A. Lang, I. G. Ryabinkin, and A. F. Izmaylov. J. Chem. Theory Comput. 2021, 17, 1, 66–78. - 3. Ç. K. Koç and S. N. Arachchige. - J. Parallel Distrib. Comput., 1991, 13, 118–122. """ import warnings +from math import acos, sin, sqrt import scipy import numpy as np +from openfermion import commutator from tangelo.toolboxes.operators.operators import QubitOperator from tangelo.toolboxes.ansatz_generator._qubit_mf import get_op_expval -def construct_acs(dis, max_ilc_gens, n_qubits): +def construct_acs(dis, n_qubits): """Driver function for constructing the anticommuting set of generators from the direct interaction set (DIS) of QCC generators. Args: dis (list of list): DIS of QCC generators. - max_ilc_gens (int): maximum number of ILC generators allowed in the ansatz. n_qubits (int): number of qubits Returns: list of QubitOperator: the anticommuting set (ACS) of ILC generators """ - bad_sln_idxs, good_sln = [], False + bad_sln_idxs = [] + n_dis_groups = len(dis) + + good_sln = False while not good_sln: - gen_idxs, ilc_gens = [idx for idx in range(max_ilc_gens) if idx not in bad_sln_idxs], [] + gen_idxs, ilc_gens = [idx for idx in range(n_dis_groups) if idx not in bad_sln_idxs], [] n_gens = len(gen_idxs) ng2, ngnq = n_gens * (n_gens + 1) // 2, n_gens * n_qubits @@ -115,9 +117,9 @@ def construct_acs(dis, max_ilc_gens, n_qubits): def gauss_elim_over_gf2(a_mat, b_vec=None): """Driver function that performs Gaussian elimination to solve A * z = b - over the binary field where b is the known solution vector. This routine - was adapted based on Ref. 3. All elements of a_mat and b_vec are assumed - to be the integers 0 or 1. + over the binary field where b is the known solution vector. The elements + of a_mat and b_vec need to be 0 or 1. If they are not provided as such, + they will be transformed accordingly. Args: a_mat (numpy array of int): rectangular matrix of dimension n x m that @@ -127,62 +129,87 @@ def gauss_elim_over_gf2(a_mat, b_vec=None): initial solution of A * z. Default, np.zeros((n, 1)). Returns: - numpy array of float: solution for the z vector of dimension (n, ) + numpy array of int: solution for the z vector of dimension (n, ) """ + # check that b_vec was provided properly; otherwise initialize as a vector of zeros n_rows, n_cols = np.shape(a_mat) - z_vals, z_sln, piv_idx = [], [-1] * n_cols, 0 - # check that b_vec was properly supplied; ortherwise initialize as a vector of zeros if not isinstance(b_vec, np.ndarray): b_vec = np.zeros((n_rows, 1)) - a_mat = np.append(a_mat, b_vec, axis=1) + + # append the initial solution vector as the last column of a_mat; update n_cols + a_mat = np.concatenate((a_mat, b_vec), axis=1).astype('int8') n_cols += 1 + + # force all entries of a_mat to be either 0 or 1. + if (abs(a_mat) > 1).any(): + warnings.warn("Reducing input matrix elements modulo 2 to create a binary matrix.", RuntimeWarning) + a_mat = (a_mat % 2).astype('int8') + + # remove duplicate rows if they exist + _, row_idxs = np.unique([tuple(row) for row in a_mat], axis=0, return_index=True) + a_mat_new = a_mat[np.sort(row_idxs)] + if a_mat_new.shape[0] != a_mat.shape[0]: + warnings.warn("Linear dependency detected in input matrix: redundant rows deleted.", RuntimeWarning) + a_mat = a_mat_new + + # remove rows of all 0s if they exist + del_idxs = [] + for i in range(a_mat.shape[0]): + if (a_mat[i][:] == 0).all(): + del_idxs.append(i) + if del_idxs: + warnings.warn("Linear dependency detected in input matrix: rows of zeros deleted.", RuntimeWarning) + a_mat = np.delete(a_mat, obj=del_idxs, axis=0) + n_rows = a_mat.shape[0] + + # begin gaussian elimination algorithm + z_sln, piv_idx = [-1]*(n_cols-1), 0 for i in range(n_cols): - a_mat_max, max_idx = 0., piv_idx - # locate the pivot index by searching each row for a non-zero value. + a_mat_max, max_idx = 0, piv_idx + # locate the pivot index by searching each row for a non-zero value for j in range(piv_idx, n_rows): - # if a pivot index is found, set the value to the col index for the row in which it was found + # if a pivot index is found, mark the value of the column index if a_mat[j, i] > a_mat_max: max_idx = j a_mat_max = a_mat[j, i] - # if a pivot index is not found in a given row, reset a_mat_max to -1 and move to the next row - elif j == n_rows-1 and a_mat_max == 0.: + # if a pivot index is not found, reset and move to the next row + elif j == n_rows-1 and a_mat_max == 0: piv_idx = max_idx - a_mat_max = -1. + a_mat_max = -1 # update the matrix by flipping the row and columns to achieve row echelon form - if a_mat_max > 0.: + if a_mat_max > 0: if max_idx > piv_idx: a_mat[[piv_idx, max_idx]] = a_mat[[max_idx, piv_idx]] for j in range(piv_idx + 1, n_rows): - if a_mat[j, i] == 1.: + if a_mat[j, i] == 1: a_mat[j, i:n_cols] = np.fmod(a_mat[j, i:n_cols] + a_mat[piv_idx, i:n_cols], 2) piv_idx += 1 - # extract the solution from the bottom to the top since it is now in row echelon form - b_vec = a_mat[0:n_rows, n_cols-1].tolist() - for i in range(n_rows - 1, -1, -1): - col_idx, z_free = -1., [] - for j in range(n_cols-1): - if a_mat[i, j] == 1.: - if col_idx == -1: - col_idx = j - else: - z_free.append(j) - if col_idx >= 0.: - z_vals.append([col_idx, z_free, b_vec[i]]) - # check for free solutions -- select 0 for the free solution - # for the ILC generator screening procedure, 0 leads to an I op and 1 leads to a Z Pauli op - for z_val in (z_vals): - b_val = z_val[2] - for z_free in (z_val[1]): - if z_sln[z_free] == -1: - z_sln[z_free] = 0. - b_val = np.fmod(b_val + z_sln[z_free], 2) - z_sln[z_val[0]] = b_val - # check that z_sln does not have any -1 values left -- if so, a solution was not found. - for z_val in z_sln: - if z_val == -1: - warnings.warn("Gaussian elimination over GF(2) failed to find a solution.", RuntimeWarning) - return np.array(z_sln) + + # extract the solution: back solve from the last row to the first + for i in range(n_rows-1, -1, -1): + # find the first non-zero coefficient for a row + j = 0 + while a_mat[i, j] == 0 and j < n_cols-2: + j += 1 + # initialize the soln then back solve + b_val = a_mat[i, n_cols-1] + if i == n_rows-1 and j == n_cols-2: + z_sln[j] = int(b_val) + else: + for k in range(j+1, n_cols-1): + if a_mat[i, k] == 1: + if z_sln[k] == -1: + z_sln[k] = 0 + else: + b_val = (b_val + z_sln[k]) % 2 + z_sln[j] = int(b_val) + + # set any remaining free variables to 0 + for i, sln_val in enumerate(z_sln): + if sln_val == -1: + z_sln[i] = 0 + return np.array(z_sln, dtype=int) def get_ilc_params_by_diag(qubit_ham, ilc_gens, qmf_var_params): @@ -215,7 +242,7 @@ def get_ilc_params_by_diag(qubit_ham, ilc_gens, qmf_var_params): # = = 1 qubit_overlap_mat[i, i] = 1. + 0j - for j in range(i + 1, n_var_params): + for j in range(i+1, n_var_params): # = = H_ji qubit_ham_mat[j, i] = get_op_expval(ilc_gens[j] * h_psi_i, qmf_var_params) @@ -234,12 +261,86 @@ def get_ilc_params_by_diag(qubit_ham, ilc_gens, qmf_var_params): gs_coefs *= -1. denom_sum, ilc_var_params = 0., [] for i in range(2): - denom_sum += pow(gs_coefs[i].real, 2.) + pow(gs_coefs[i].imag, 2.) + denom_sum += abs(gs_coefs[i])**2 beta_1 = np.arcsin(gs_coefs[1] / np.sqrt(denom_sum)) ilc_var_params.append(beta_1.real) for i in range(2, n_var_params): - denom_sum += pow(gs_coefs[i].real, 2.) + pow(gs_coefs[i].imag, 2.) + denom_sum += abs(gs_coefs[i])**2 beta = np.arcsin(gs_coefs[i] / np.sqrt(denom_sum)) ilc_var_params.append(beta.real) del ilc_gens[0] return ilc_var_params + + +def build_ilc_qubit_op_list(acs_gens, ilc_params): + """Returns a list of 2N - 1 ILC generators to facilitate generation of a circuit + based on symmetric Trotter-Suzuki decomposition. The ILC generators are ordered + according to Eq. C1 in Appendix C of Ref. 1. + + Args: + acs_gens (list of QubitOperator): The list of ILC Pauli word generators + selected from characteristic ACS groups. + ilc_params (list or numpy array of float): The ILC variational parameters + arranged such that their ordering matches the order of acs_gens. + + Returns: + list of QubitOperator: list of ILC ansatz operator generators. + """ + + n_amps = len(ilc_params) + ilc_op_list = [-.5 * ilc_params[i] * acs_gens[i] for i in range(n_amps-1, 0, -1)] + ilc_op_list += [ilc_params[0] * acs_gens[0]] + ilc_op_list += [-.5 * ilc_params[i] * acs_gens[i] for i in range(1, n_amps)] + return ilc_op_list + + +def ilc_op_dress(qubit_op, ilc_gens, ilc_params): + """Performs transformation of a qubit operator with the ACS of ILC generators and + parameters. For a set of N generators, each qubit operator transformation results + in quadratic (N * (N-1) / 2) growth of the number of its terms. + + Args: + qubit_op (QubitOperator): A qubit operator to be dressed. + ilc_gens (list of QubitOperator): The list of ILC Pauli word generators + selected from a user-specified number of characteristic ACS groups. + ilc_params (list or numpy array of float): The ILC variational parameters + arranged such that their ordering matches the ordering of ilc_gens. + + Returns: + QubitOperator: Dressed qubit operator. + """ + + # first, recast the beta parameters into the set of coefficients {c_n} + n_amps = len(ilc_params) + coef_norm = 1. + coefs = [0.] * n_amps + + # See Ref. 1, Appendix C, Eqs. C3 and C4: + # sin b_n = c_n; sin_b_n-1 = c_n-1 / sqrt(1-|c_n|**2); + # sin_b_n-2 = c_n-2 / sqrt(1-|c_n|**2-|c_n-1|**2) ... + for i in range(n_amps-1, -1, -1): + coef = sqrt(coef_norm) * sin(ilc_params[i]) + coefs[i] = coef + coef_norm -= coef**2 + + # the remainder of coef_norm is |c_0|^2 + coefs.insert(0, -sqrt(coef_norm)) + + # second, recast {c_n} into tau, {alpha_i}; + # c_0 = cos(tau); c_n = sin(tau) * alpha_n for n > 0 + tau = acos(coefs[0]) + sin_tau = sin(tau) + alphas = [coefs[i]/sin_tau for i in range(1, n_amps+1)] + + # third, dress the qubit operator according to Eqs. 17, 18 in Ref. 2 + sin2_tau = sin_tau**2 + sin_2tau = sin(2.*tau) + qop_dress = coefs[0]**2 * qubit_op + for i in range(n_amps): + qop_dress += sin2_tau * alphas[i]**2 * ilc_gens[i] * qubit_op * ilc_gens[i]\ + - .5j * sin_2tau * alphas[i] * commutator(qubit_op, ilc_gens[i]) + for j in range(i+1, n_amps): + qop_dress += sin2_tau * alphas[i] * alphas[j] * (ilc_gens[i] * qubit_op * ilc_gens[j] + + ilc_gens[j] * qubit_op * ilc_gens[i]) + qop_dress.compress() + return qop_dress diff --git a/tangelo/toolboxes/ansatz_generator/ilc.py b/tangelo/toolboxes/ansatz_generator/ilc.py index ba708a0ef..167310df4 100755 --- a/tangelo/toolboxes/ansatz_generator/ilc.py +++ b/tangelo/toolboxes/ansatz_generator/ilc.py @@ -34,11 +34,13 @@ from tangelo.toolboxes.qubit_mappings.mapping_transform import get_qubit_number,\ fermion_to_qubit_mapping from tangelo.linq import Circuit +from tangelo import SecondQuantizedMolecule from tangelo.toolboxes.ansatz_generator.ansatz import Ansatz from tangelo.toolboxes.ansatz_generator.ansatz_utils import exp_pauliword_to_gates from tangelo.toolboxes.ansatz_generator._qubit_mf import init_qmf_from_hf, get_qmf_circuit, purify_qmf_state from tangelo.toolboxes.ansatz_generator._qubit_cc import construct_dis -from tangelo.toolboxes.ansatz_generator._qubit_ilc import construct_acs, get_ilc_params_by_diag +from tangelo.toolboxes.ansatz_generator._qubit_ilc import construct_acs, get_ilc_params_by_diag,\ + build_ilc_qubit_op_list class ILC(Ansatz): @@ -49,7 +51,9 @@ class ILC(Ansatz): state is obtained using a RHF or ROHF Hamiltonian, respectively. Args: - molecule (SecondQuantizedMolecule): The molecular system. + molecule (SecondQuantizedMolecule or dict): The molecular system, which can + be passed as a SecondQuantizedMolecule or a dictionary with keys that + specify n_spinoribtals, n_electrons, and spin. Default, None. mapping (str): One of the supported mapping identifiers. Default, "jw". up_then_down (bool): Change basis ordering putting all spin-up orbitals first, followed by all spin-down. Default, False. @@ -75,10 +79,19 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, acs=None, qmf_circuit=None, qmf_var_params=None, qubit_ham=None, ilc_tau_guess=1e-2, deilc_dtau_thresh=1e-3, max_ilc_gens=None): - if not molecule: - raise ValueError("An instance of SecondQuantizedMolecule is required for initializing " - "the self.__class__.__name__ ansatz class.") + if not molecule and not (isinstance(molecule, SecondQuantizedMolecule) and isinstance(molecule, dict)): + raise ValueError("An instance of SecondQuantizedMolecule or a dict is required for " + "initializing the self.__class__.__name__ ansatz class.") self.molecule = molecule + if isinstance(self.molecule, SecondQuantizedMolecule): + self.n_spinorbitals = self.molecule.n_active_sos + self.n_electrons = self.molecule.n_electrons + self.spin = self.molecule.spin + elif isinstance(self.molecule, dict): + self.n_spinorbitals = self.molecule["n_spinorbitals"] + self.n_electrons = self.molecule["n_electrons"] + self.spin = self.molecule["spin"] + self.mapping = mapping self.up_then_down = up_then_down if self.mapping.lower() == "jw" and not self.up_then_down: @@ -87,31 +100,24 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, acs=None, "with the self.__class__.__name__ ansatz.", RuntimeWarning) self.up_then_down = True - self.n_spinorbitals = self.molecule.n_active_sos if self.n_spinorbitals % 2 != 0: raise ValueError("The total number of spin-orbitals should be even.") - - self.spin = molecule.spin - self.fermi_ham = self.molecule.fermionic_hamiltonian - self.n_electrons = self.molecule.n_electrons self.n_qubits = get_qubit_number(self.mapping, self.n_spinorbitals) self.qubit_ham = qubit_ham - if qubit_ham is None: + if not qubit_ham: self.fermi_ham = self.molecule.fermionic_hamiltonian self.qubit_ham = fermion_to_qubit_mapping(self.fermi_ham, self.mapping, self.n_spinorbitals, self.n_electrons, self.up_then_down, self.spin) - self.qmf_var_params = qmf_var_params - if self.qmf_var_params is None: + self.qmf_var_params = np.array(qmf_var_params) if isinstance(qmf_var_params, list) else qmf_var_params + if not isinstance(self.qmf_var_params, np.ndarray): self.qmf_var_params = init_qmf_from_hf(self.n_spinorbitals, self.n_electrons, self.mapping, self.up_then_down, self.spin) - elif isinstance(self.qmf_var_params, list): - self.qmf_var_params = np.array(self.qmf_var_params) if self.qmf_var_params.size != 2 * self.n_qubits: raise ValueError("The number of QMF variational parameters must be 2 * n_qubits.") - + self.n_qmf_params = 2 * self.n_qubits self.qmf_circuit = qmf_circuit self.acs = acs @@ -119,18 +125,10 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, acs=None, self.deilc_dtau_thresh = deilc_dtau_thresh self.max_ilc_gens = max_ilc_gens - # Get purified QMF parameters and build the DIS & ACS or use a list of generators. - if self.acs is None: - pure_var_params = purify_qmf_state(self.qmf_var_params, self.n_spinorbitals, - self.n_electrons, self.mapping, self.up_then_down, self.spin) - self.dis = construct_dis(self.qubit_ham, pure_var_params, self.deilc_dtau_thresh) - self.max_ilc_gens = len(self.dis) if self.max_ilc_gens is None\ - else min(len(self.dis), self.max_ilc_gens) - self.acs = construct_acs(self.dis, self.max_ilc_gens, self.n_qubits) - self.n_var_params = len(self.acs) - else: - self.dis = None - self.n_var_params = len(self.acs) + # Build the ACS of ILC generators + self.n_ilc_params = 0 + self._get_ilc_generators() + self.n_var_params = self.n_qmf_params + self.n_ilc_params # Supported reference state initialization self.supported_reference_state = {"HF"} @@ -142,7 +140,6 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, acs=None, self.var_params_default = "diag" self.var_params = None self.rebuild_dis = False - self.rebuild_acs = False self.ilc_circuit = None self.circuit = None @@ -172,30 +169,37 @@ def set_var_params(self, var_params=None): # Initialize ILC parameters by matrix diagonalization (see Appendix B, Refs. 1 & 2). elif var_params == "diag": initial_var_params = get_ilc_params_by_diag(self.qubit_ham, self.acs, self.qmf_var_params) + # Insert the QMF variational parameters at the beginning. + initial_var_params = np.concatenate((self.qmf_var_params, initial_var_params)) else: initial_var_params = np.array(var_params) - if initial_var_params.size != self.n_var_params: - raise ValueError(f"Expected {self.n_var_params} variational parameters but " - f"received {initial_var_params.size}.") self.var_params = initial_var_params + if initial_var_params.size != self.n_var_params: + raise ValueError(f"Expected {self.n_var_params} variational parameters but " + f"received {initial_var_params.size}.") return initial_var_params 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 operator. """ + with the transform used to obtain the qubit operator. """ if self.default_reference_state not in self.supported_reference_state: raise ValueError(f"Only supported reference state methods are: " f"{self.supported_reference_state}.") if self.default_reference_state == "HF": - reference_state_circuit = get_qmf_circuit(self.qmf_var_params, False) + reference_state_circuit = get_qmf_circuit(self.qmf_var_params, True) return reference_state_circuit def build_circuit(self, var_params=None): """Build and return the quantum circuit implementing the state preparation ansatz (with currently specified initial_state and var_params). """ + # Build the DIS or specify a list of generators; updates the number of QCC parameters + self._get_ilc_generators() + self.n_var_params = self.n_qmf_params + self.n_ilc_params + + # Get the variational parameters needed for the QCC unitary operator and circuit if var_params is not None: self.set_var_params(var_params) elif self.var_params is None: @@ -205,12 +209,12 @@ def build_circuit(self, var_params=None): if self.qmf_circuit is None: self.qmf_circuit = self.prepare_reference_state() - # Build create the list of ILC qubit operators - self.ilc_op_list = self._get_ilc_op() + # Build the QCC ansatz qubit operator + ilc_op_list = build_ilc_qubit_op_list(self.acs, self.var_params[2*self.n_qubits:]) - # Obtain quantum circuit through trotterization of the list of ILC operators + # Obtain quantum circuit corresponding to the list of ILC generators pauli_word_gates = [] - for ilc_op in self.ilc_op_list: + for ilc_op in ilc_op_list: pauli_word, coef = list(ilc_op.terms.items())[0] pauli_word_gates += exp_pauliword_to_gates(pauli_word, coef, variational=True) self.ilc_circuit = Circuit(pauli_word_gates) @@ -226,39 +230,29 @@ def update_var_params(self, var_params): self.set_var_params(var_params) # Build the ILC ansatz operator - self.ilc_op_list = self._get_ilc_op() + ilc_op_list = build_ilc_qubit_op_list(self.acs, self.var_params[2*self.n_qubits:]) pauli_word_gates = [] - for ilc_op in self.ilc_op_list: + for ilc_op in ilc_op_list: pauli_word, coef = list(ilc_op.terms.items())[0] pauli_word_gates += exp_pauliword_to_gates(pauli_word, coef, variational=True) self.ilc_circuit = Circuit(pauli_word_gates) self.circuit = self.qmf_circuit + self.ilc_circuit if self.qmf_circuit.size != 0\ else self.ilc_circuit - def _get_ilc_op(self): - """Returns the ILC operators ordered according to the argument of - Eq. C1, Appendix C, Ref. 1. - - Returns: - list of QubitOperator: the list of ILC qubit operators - """ + def _get_ilc_generators(self): + """ Prepares the ILC ansatz by purifying the QMF state, constructing the ACS, + and selecting representative generators from the top candidate ACS groups. """ - # Rebuild DIS & ACS in case qubit_ham changed or they and qubit_op_list don't exist - if self.rebuild_dis or self.rebuild_acs or not self.acs: - pure_var_params = purify_qmf_state(self.qmf_var_params, self.n_spinorbitals, - self.n_electrons, self.mapping, self.up_then_down, self.spin) + if not self.acs: + pure_var_params = purify_qmf_state(self.qmf_var_params, self.n_spinorbitals, self.n_electrons, + self.mapping, self.up_then_down, self.spin) self.dis = construct_dis(self.qubit_ham, pure_var_params, self.deilc_dtau_thresh) - self.max_ilc_gens = len(self.dis) if self.max_ilc_gens is None\ - else min(len(self.dis), self.max_ilc_gens) - self.acs = construct_acs(self.dis, self.max_ilc_gens, self.n_qubits) - self.ilc_op_list = None - - # Build the ILC qubit operator list - ilc_op_list = [] - for i in range(self.n_var_params - 1, 0, -1): - ilc_op_list.append(-0.5 * self.var_params[i] * self.acs[i]) - ilc_op_list.append(-self.var_params[0] * self.acs[0]) - for i in range(1, self.n_var_params): - ilc_op_list.append(-0.5 * self.var_params[i] * self.acs[i]) - return ilc_op_list + self.acs = construct_acs(self.dis, self.n_qubits) + if self.max_ilc_gens: + self.n_ilc_params = min(len(self.acs), self.max_ilc_gens) + del self.acs[self.n_ilc_params:] + else: + self.n_ilc_params = len(self.acs) + else: + self.n_ilc_params = len(self.acs) diff --git a/tangelo/toolboxes/ansatz_generator/qcc.py b/tangelo/toolboxes/ansatz_generator/qcc.py index 6ee7b9863..3801061ba 100755 --- a/tangelo/toolboxes/ansatz_generator/qcc.py +++ b/tangelo/toolboxes/ansatz_generator/qcc.py @@ -40,6 +40,7 @@ from tangelo.toolboxes.qubit_mappings.mapping_transform import get_qubit_number,\ fermion_to_qubit_mapping from tangelo.linq import Circuit +from tangelo import SecondQuantizedMolecule from tangelo.toolboxes.ansatz_generator.ansatz import Ansatz from tangelo.toolboxes.ansatz_generator.ansatz_utils import exp_pauliword_to_gates from tangelo.toolboxes.ansatz_generator._qubit_mf import init_qmf_from_hf, get_qmf_circuit, purify_qmf_state @@ -54,7 +55,9 @@ class QCC(Ansatz): state is obtained using a RHF or ROHF Hamiltonian, respectively. Args: - molecule (SecondQuantizedMolecule): The molecular system. + molecule (SecondQuantizedMolecule or dict): The molecular system, which can + be passed as a SecondQuantizedMolecule or a dictionary with keys that + specify n_spinoribtals, n_electrons, and spin. Default, None. mapping (str): One of the supported qubit mapping identifiers. Default, "jw". up_then_down (bool): Change basis ordering putting all spin-up orbitals first, followed by all spin-down. Default, False. @@ -80,10 +83,19 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, dis=None, qmf_circuit=None, qmf_var_params=None, qubit_ham=None, qcc_tau_guess=1e-2, deqcc_dtau_thresh=1e-3, max_qcc_gens=None): - if not molecule: - raise ValueError("An instance of SecondQuantizedMolecule is required for initializing " - "the self.__class__.__name__ ansatz class.") + if not molecule and not (isinstance(molecule, SecondQuantizedMolecule) and isinstance(molecule, dict)): + raise ValueError("An instance of SecondQuantizedMolecule or a dict is required for " + "initializing the self.__class__.__name__ ansatz class.") self.molecule = molecule + if isinstance(self.molecule, SecondQuantizedMolecule): + self.n_spinorbitals = self.molecule.n_active_sos + self.n_electrons = self.molecule.n_electrons + self.spin = self.molecule.spin + elif isinstance(self.molecule, dict): + self.n_spinorbitals = self.molecule["n_spinorbitals"] + self.n_electrons = self.molecule["n_electrons"] + self.spin = self.molecule["spin"] + self.mapping = mapping self.up_then_down = up_then_down if self.mapping.lower() == "jw" and not self.up_then_down: @@ -92,13 +104,8 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, dis=None, "with the self.__class__.__name__ ansatz.", RuntimeWarning) self.up_then_down = True - self.n_spinorbitals = self.molecule.n_active_sos if self.n_spinorbitals % 2 != 0: raise ValueError("The total number of spin-orbitals should be even.") - - self.spin = molecule.spin - self.fermi_ham = self.molecule.fermionic_hamiltonian - self.n_electrons = self.molecule.n_electrons self.n_qubits = get_qubit_number(self.mapping, self.n_spinorbitals) self.qubit_ham = qubit_ham @@ -108,12 +115,10 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, dis=None, self.n_spinorbitals, self.n_electrons, self.up_then_down, self.spin) - self.qmf_var_params = qmf_var_params - if not self.qmf_var_params: + self.qmf_var_params = np.array(qmf_var_params) if isinstance(qmf_var_params, list) else qmf_var_params + if not isinstance(self.qmf_var_params, np.ndarray): self.qmf_var_params = init_qmf_from_hf(self.n_spinorbitals, self.n_electrons, self.mapping, self.up_then_down, self.spin) - elif isinstance(self.qmf_var_params, list): - self.qmf_var_params = np.array(self.qmf_var_params) if self.qmf_var_params.size != 2 * self.n_qubits: raise ValueError("The number of QMF variational parameters must be 2 * n_qubits.") self.n_qmf_params = 2 * self.n_qubits @@ -125,6 +130,7 @@ def __init__(self, molecule, mapping="jw", up_then_down=False, dis=None, self.max_qcc_gens = max_qcc_gens # Build the DIS or specify a list of generators; updates the number of QCC parameters + self.n_qcc_params = 0 self._get_qcc_generators() self.n_var_params = self.n_qmf_params + self.n_qcc_params diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py b/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py index 706862ba4..334295067 100644 --- a/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py +++ b/tangelo/toolboxes/ansatz_generator/tests/test_ilc.py @@ -39,21 +39,18 @@ def test_ilc_set_var_params(): ilc_ansatz = ILC(mol_H2_sto3g, up_then_down=True) - one_zero = np.zeros((1,), dtype=float) + nine_zeros = np.zeros((9,), dtype=float) - ilc_ansatz.set_var_params("qmf_state") - np.testing.assert_array_almost_equal(ilc_ansatz.var_params, one_zero, decimal=6) + ilc_ansatz.set_var_params([0.] * 9) + np.testing.assert_array_almost_equal(ilc_ansatz.var_params, nine_zeros, decimal=6) - ilc_ansatz.set_var_params([0.]) - np.testing.assert_array_almost_equal(ilc_ansatz.var_params, one_zero, decimal=6) + nine_tenths = 0.1 * np.ones((9,)) - one_tenth = 0.1 * np.ones((1,)) + ilc_ansatz.set_var_params([0.1] * 9) + np.testing.assert_array_almost_equal(ilc_ansatz.var_params, nine_tenths, decimal=6) - ilc_ansatz.set_var_params([0.1]) - np.testing.assert_array_almost_equal(ilc_ansatz.var_params, one_tenth, decimal=6) - - ilc_ansatz.set_var_params(np.array([0.1])) - np.testing.assert_array_almost_equal(ilc_ansatz.var_params, one_tenth, decimal=6) + ilc_ansatz.set_var_params(np.array([0.1] * 9)) + np.testing.assert_array_almost_equal(ilc_ansatz.var_params, nine_tenths, decimal=6) def test_ilc_incorrect_number_var_params(self): """ Return an error if user provide incorrect number of variational parameters """ @@ -64,7 +61,7 @@ def test_ilc_incorrect_number_var_params(self): @staticmethod def test_gauss_elim_over_gf2_sqrmat(): - """ Verify behavior of the Gaussian elimination over the binary field function. """ + """ Verify behavior of the Gaussian elimination for a square matrix.""" # a_matrix stores the action of A * z over GF(2); dimension is n x m a_matrix = np.array([[1, 0, 1, 1], [1, 1, 0, 1], [1, 1, 1, 0], [0, 0, 1, 0]]) @@ -82,7 +79,7 @@ def test_gauss_elim_over_gf2_sqrmat(): @staticmethod def test_gauss_elim_over_gf2_rectmat(): - """ Verify behavior of the Gaussian elimination over the binary field function. """ + """ Verify behavior of the Gaussian elimination for a rectangular matrix.""" # a_matrix stores the action of A * z over GF(2); dimension is n x m a_matrix = np.array([[0, 0, 1, 0, 1], [1, 1, 0, 0, 0], [0, 0, 0, 1, 1]]) @@ -100,16 +97,36 @@ def test_gauss_elim_over_gf2_rectmat(): @staticmethod def test_gauss_elim_over_gf2_lindep(): - """ Verify behavior of the Gaussian elimination over the binary field function. """ + """ Verify behavior of the Gaussian elimination when linear dependence + in the form of duplicate rows arises in a matrix.""" # a_matrix stores the action of A * z over GF(2); dimension is n x m a_matrix = np.array([[0, 0, 1, 0, 1], [0, 0, 1, 0, 1], [0, 0, 0, 1, 1]]) # b_vec stores the solution vector for the equation A * z = b_vec; dimension is n x 1 - b_vec = np.array([1, 0, 1]).reshape((3, 1)) + b_vec = np.array([0, 0, 1]).reshape((3, 1)) + + # z_ref stores the serves as the reference for the output of gauss_elim_over_gf2 + z_ref = np.array([0, 0, 0, 1, 0]) + + # solve A * z = b and compare to reference solution + z_sln = gauss_elim_over_gf2(a_matrix, b_vec) + + np.testing.assert_array_almost_equal(z_sln, z_ref, decimal=6) + + @staticmethod + def test_gauss_elim_over_gf2_lindep2(): + """ Verify behavior of the Gaussian elimination when linear dependence + in the form of a row of zeros arises in a matrix.""" + + # a_matrix stores the action of A * z over GF(2); dimension is n x m + a_matrix = np.array([[0, 1, 0, 1, 0], [1, 0, 1, 0, 1], [0, 0, 0, 0, 0]]) + + # b_vec stores the solution vector for the equation A * z = b_vec; dimension is n x 1 + b_vec = np.array([1, 1, 0]).reshape((3, 1)) # z_ref stores the serves as the reference for the output of gauss_elim_over_gf2 - z_ref = np.array([-1, -1, 1, 1, 0]) + z_ref = np.array([1, 1, 0, 0, 0]) # solve A * z = b and compare to reference solution z_sln = gauss_elim_over_gf2(a_matrix, b_vec) @@ -117,44 +134,52 @@ def test_gauss_elim_over_gf2_lindep(): np.testing.assert_array_almost_equal(z_sln, z_ref, decimal=6) def test_ilc_h2(self): - """ Verify closed-shell functionality when using the ILC class separately for H2 """ + """ Verify closed-shell functionality when using the ILC class separately for H2.""" - # Specify the mutually anticommuting set (ACS) of ILC generators and parameters. - acs = [QubitOperator("X0 Y1 Y2 Y3")] - # The QMF parameters are automatically set if the argument qmf_var_params is not given. - ilc_var_params = [0.11360304] - ilc_ansatz = ILC(mol_H2_sto3g, "jw", True, acs) + # Specify the qubit operators from the anticommuting set (ACS) of ILC generators. + acs = [QubitOperator("Y0 X1")] + ilc_ansatz = ILC(mol_H2_sto3g, mapping="scbk", up_then_down=True, acs=acs) - # Build the combined QMF (determined automatically if not specified) + ILC circuit. + # Build the ILC circuit, which is prepended by the qubit mean field (QMF) circuit. ilc_ansatz.build_circuit() # Get qubit hamiltonian for energy evaluation qubit_hamiltonian = ilc_ansatz.qubit_ham + # The QMF and ILC parameters can both be specified; determined automatically otherwise. + qmf_var_params = [3.14159265e+00, 3.14159265e+00, -7.59061327e-12, 0.] + ilc_var_params = [1.12894599e-01] + var_params = qmf_var_params + ilc_var_params + ilc_ansatz.update_var_params(var_params) # Assert energy returned is as expected for given parameters - ilc_ansatz.update_var_params(ilc_var_params) energy = sim.get_expectation_value(qubit_hamiltonian, ilc_ansatz.circuit) - self.assertAlmostEqual(energy, -1.1372697, delta=1e-6) + self.assertAlmostEqual(energy, -1.137270126, delta=1e-6) def test_ilc_h4_cation(self): """ Verify restricted open-shell functionality when using the ILC class for H4+ """ - # Specify the mutually anticommuting set (ACS) of ILC generators and parameters. - acs = [QubitOperator("Y0 Z2 X4 Z6"), QubitOperator("Y1 Y2 Z4 X5 Y6"), QubitOperator("X0 Z2 Z4 Y6"), - QubitOperator("X1 Y2 X4 Z6"), QubitOperator("Y1 Y2 X4 Y5 Z6"), QubitOperator("Y1 Y2 Z4 Z5 Y6"), - QubitOperator("Y0 Z1 Z2 Y5 Y6"), QubitOperator("Y0 Z1 Z2 Y4 Y5 Z6")] - # The QMF parameters are automatically set if the argument qmf_var_params is not given. - ilc_var_params = [ 0.14017492, -0.10792805, -0.05835484, 0.12468933, 0.07173118, 0.04683807, 0.02852163, -0.03133538] - ilc_ansatz = ILC(mol_H4_cation_sto3g, "bk", False, acs) + # Specify the qubit operators from the anticommuting set (ACS) of ILC generators. + acs = [QubitOperator("X0 X1 X2 X3 X4 Y5"), QubitOperator("X1 Z2 X3 X4 Y5"), + QubitOperator("Y0 X1 Z2 X3 X4 Z5"), QubitOperator("Z0 X1 X2 X3 X4 Y5"), + QubitOperator("X1 Y2 X3 X4"), QubitOperator("Y1 X3 X4"), + QubitOperator("Y0 X1 Z2 X3 X4 X5"), QubitOperator("Y0 X1 X2 X3 X4")] + ilc_ansatz = ILC(mol_H4_cation_sto3g, mapping="scbk", up_then_down=True, acs=acs) - # Build the combined QMF (determined automatically if not specified) + ILC circuit. + # Build the ILC circuit, which is prepended by the qubit mean field (QMF) circuit. ilc_ansatz.build_circuit() # Get qubit hamiltonian for energy evaluation qubit_hamiltonian = ilc_ansatz.qubit_ham + # The QMF and ILC parameters can both be specified; determined automatically otherwise. + qmf_var_params = [ 3.14159265e+00, -1.02576971e-11, 1.35522331e-11, 3.14159265e+00, + 3.14159265e+00, -5.62116001e-11, -1.41419277e-11, -2.36789365e-11, + -5.53225030e-11, -3.56400157e-11, -2.61030058e-11, -3.55652002e-11] + ilc_var_params = [ 0.14001419, -0.10827113, 0.05840200, -0.12364925, + -0.07275071, -0.04703495, 0.02925292, 0.03145765] + var_params = qmf_var_params + ilc_var_params # Assert energy returned is as expected for given parameters - ilc_ansatz.update_var_params(ilc_var_params) + ilc_ansatz.update_var_params(var_params) energy = sim.get_expectation_value(qubit_hamiltonian, ilc_ansatz.circuit) self.assertAlmostEqual(energy, -1.6379638, delta=1e-6) diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py b/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py index 0327f3ff1..d175b643b 100644 --- a/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py +++ b/tangelo/toolboxes/ansatz_generator/tests/test_qcc.py @@ -74,7 +74,7 @@ def test_qcc_h2(self): # Get qubit hamiltonian for energy evaluation qubit_hamiltonian = qcc_ansatz.qubit_ham - # The QMF and QCC parameters can both be specified; determined automatically othersise. + # The QMF and QCC parameters can both be specified; determined automatically otherwise. qmf_var_params = [ 3.14159265e+00, -2.42743256e-08, 3.14159266e+00, -3.27162543e-08, 3.08514545e-09, 3.08514545e-09, 3.08514545e-09, 3.08514545e-09] qcc_var_params = [-2.26136280e-01] @@ -101,7 +101,7 @@ def test_qmf_qcc_h4_cation(self): # Get qubit hamiltonian for energy evaluation qubit_hamiltonian = qcc_ansatz.qubit_ham - # The QMF and QCC parameters can both be specified; determined automatically othersise. + # The QMF and QCC parameters can both be specified; determined automatically otherwise. qmf_var_params = [3.14159302e+00, 6.20193478e-07, 1.51226426e-06, 3.14159350e+00, 3.14159349e+00, 7.88310582e-07, 3.96032530e+00, 2.26734374e+00, 3.22127001e+00, 5.77997401e-01, 5.51422406e+00, 6.26513711e+00] @@ -128,7 +128,7 @@ def test_qmf_qcc_h4_double_cation(self): # Get qubit hamiltonian for energy evaluation qubit_hamiltonian = load_operator("mol_H4_doublecation_minao_qubitham_bk_updown.data", data_directory=pwd_this_test+"/data", plain_text=True) - # The QMF and QCC parameters can both be specified; determined automatically othersise. + # The QMF and QCC parameters can both be specified; determined automatically otherwise. qmf_var_params = [3.14159247e+00, 3.14158884e+00, 1.37660700e-06, 3.14159264e+00, 3.14159219e+00, 3.14158908e+00, 0.00000000e+00, 0.00000000e+00, 6.94108155e-01, 1.03928030e-01, 5.14029803e+00, 2.81850365e+00, diff --git a/tangelo/toolboxes/operators/operators.py b/tangelo/toolboxes/operators/operators.py index a1a8330cf..7b28f21f9 100644 --- a/tangelo/toolboxes/operators/operators.py +++ b/tangelo/toolboxes/operators/operators.py @@ -19,6 +19,8 @@ from math import sqrt from collections import OrderedDict +from scipy.special import comb + # Later on, if needed, we can extract the code for the operators themselves to remove the dependencies and customize import openfermion @@ -65,6 +67,21 @@ def frobenius_norm_compression(self, epsilon, n_qubits): self.terms = compressed_op self.compress() + def get_max_number_hamiltonian_terms(self, n_qubits): + """Compute the possible number of terms for a qubit Hamiltonian. In the + absence of an external magnetic field, each Hamiltonian term must have + an even number of Pauli Y operators to preserve time-reversal symmetry. + See J. Chem. Theory Comput. 2020, 16, 2, 1055–1063 for more details. + + Args: + n_qubits (int): Number of qubits in the register. + + Returns: + int: The maximum number of possible qubit Hamiltonian terms. + """ + + return sum([comb(n_qubits, 2*i, exact=True) * 3**(n_qubits-2*i) for i in range(n_qubits//2)]) + class QubitHamiltonian(QubitOperator): """QubitHamiltonian objects are essentially openfermion.QubitOperator