diff --git a/.github/workflows/ci_backends.yml b/.github/workflows/ci_backends.yml index c9c5852c..9ae5d9ab 100644 --- a/.github/workflows/ci_backends.yml +++ b/.github/workflows/ci_backends.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7, 3.8] + python-version: ['3.9'] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/ci_basic.yml b/.github/workflows/ci_basic.yml index d57d144b..80551785 100644 --- a/.github/workflows/ci_basic.yml +++ b/.github/workflows/ci_basic.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7, 3.8] + python-version: ['3.9', '3.10'] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/ci_basic_autograd.yml b/.github/workflows/ci_basic_autograd.yml index bfcb7dfc..25799934 100644 --- a/.github/workflows/ci_basic_autograd.yml +++ b/.github/workflows/ci_basic_autograd.yml @@ -1,7 +1,7 @@ # This workflow will install Python dependencies, run tests and lint with a variety of Python versions # For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions -name: Tequila-Test-Basic-Autograd +name: Tequila-Test-Basic-JAX on: push: @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7, 3.8, 3.9] + python-version: ['3.9', '3.10'] steps: - uses: actions/checkout@v2 @@ -28,8 +28,8 @@ jobs: run: | python -m pip install --upgrade pip pip install -r requirements.txt - pip uninstall -y jax jaxlib - pip install autograd + pip uninstall autograd -y + pip install jax jaxlib - name: Lint with flake8 run: | pip install flake8 diff --git a/.github/workflows/ci_chemistry_madness.yml b/.github/workflows/ci_chemistry_madness.yml index 2dd7e736..7c903227 100644 --- a/.github/workflows/ci_chemistry_madness.yml +++ b/.github/workflows/ci_chemistry_madness.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7] + python-version: ['3.10'] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/ci_chemistry_psi4.yml b/.github/workflows/ci_chemistry_psi4.yml index ef9c1e3a..98481972 100644 --- a/.github/workflows/ci_chemistry_psi4.yml +++ b/.github/workflows/ci_chemistry_psi4.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.8] + python-version: [3.9] steps: - uses: actions/checkout@v2 @@ -50,7 +50,7 @@ jobs: source $HOME/.bashrc source $CONDABASE/bin/activate conda activate test_psi4 - conda install psi4 python=3.8 -c psi4 + conda install psi4 python=3.9 -c conda-forge python -m pip install --upgrade pip python -m pip install -r requirements.txt python -m pip install -e . diff --git a/.github/workflows/ci_chemistry_pyscf.yml b/.github/workflows/ci_chemistry_pyscf.yml index 64d75b44..f68e6c64 100644 --- a/.github/workflows/ci_chemistry_pyscf.yml +++ b/.github/workflows/ci_chemistry_pyscf.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7] + python-version: [3.9] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/ci_conda_madness.yml b/.github/workflows/ci_conda_madness.yml index fc70be6d..ed042d4a 100644 --- a/.github/workflows/ci_conda_madness.yml +++ b/.github/workflows/ci_conda_madness.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7] + python-version: [3.9] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/ci_ml.yml b/.github/workflows/ci_ml.yml index aeb00d54..e649b377 100644 --- a/.github/workflows/ci_ml.yml +++ b/.github/workflows/ci_ml.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7, 3.8] + python-version: [3.9] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/ci_optimizers.yml b/.github/workflows/ci_optimizers.yml deleted file mode 100644 index 0d544f07..00000000 --- a/.github/workflows/ci_optimizers.yml +++ /dev/null @@ -1,40 +0,0 @@ -# This workflow will install Python dependencies, run tests and lint with a variety of Python versions -# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions - -name: Tequila-Test-Optimizers - -on: - push: - branches: [ master, devel ] - pull_request: - branches: [ master, devel ] - -jobs: - - build: - - runs-on: ubuntu-latest - strategy: - matrix: - python-version: [3.7, 3.8] - - steps: - - uses: actions/checkout@v2 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v1 - with: - python-version: ${{ matrix.python-version }} - - name: Install quantum backends - run: | - pip install "cirq" "qiskit>=0.30" "qulacs" - pip uninstall pyquil -y - - - name: Install and test GPyOpt interface (no qibo) - run: | - python -m pip install --upgrade pip - pip install --upgrade -r requirements.txt - pip install --upgrade -r requirements_gpyopt.txt - pip install -e . - - pytest tests/test_gpyopt.py --slow - diff --git a/.github/workflows/ci_pyquil.yml b/.github/workflows/ci_pyquil.yml index 92e39757..926ec8dc 100644 --- a/.github/workflows/ci_pyquil.yml +++ b/.github/workflows/ci_pyquil.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.7] + python-version: [3.9] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml index 4f3f6f3e..09e44946 100644 --- a/.github/workflows/pypi.yml +++ b/.github/workflows/pypi.yml @@ -30,7 +30,7 @@ jobs: # remove qulacs from dependencies (issues with windows and mac) # users need to install themselves if they want it - cat requirements.txt | sed "s|qulacs|#qulacs|g" > tmp.txt + #cat requirements.txt | sed "s|qulacs|#qulacs|g" > tmp.txt rm requirements.txt mv tmp.txt requirements.txt python setup.py sdist bdist_wheel diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml deleted file mode 100644 index 1c80674b..00000000 --- a/.github/workflows/release.yml +++ /dev/null @@ -1,22 +0,0 @@ -name: 'Tweet when released' -on: - release: - types: [released] - -jobs: - tweet: - runs-on: ubuntu-latest - steps: - - name: Tweet - id: tweet - uses: doomspec/auto-tweet-v2@v0.1.0 - env: - CONSUMER_API_KEY: ${{ secrets.TWITTER_CONSUMER_API_KEY }} - CONSUMER_API_SECRET_KEY: ${{ secrets.TWITTER_CONSUMER_API_SECRET_KEY }} - ACCESS_TOKEN: ${{ secrets.TWITTER_ACCESS_TOKEN }} - ACCESS_TOKEN_SECRET: ${{ secrets.TWITTER_ACCESS_TOKEN_SECRET }} - with: - text: | - New version released: ${{ github.event.release.name }} - ${{ github.event.release.html_url }} - - run: echo ${{ steps.tweet.outputs.response }} diff --git a/README.md b/README.md index 3cec3281..63c528f7 100644 --- a/README.md +++ b/README.md @@ -14,7 +14,7 @@ Tequila can execute the underlying quantum expectation values on state of the ar - [talks and slides](https://kottmanj.github.io/talks_and_material/) # Installation -Recommended Python version is 3.8-3.9. +Recommended Python version is 3.9 (3.10). Tequila supports linux, osx and windows. However, not all optional dependencies are supported on windows. ## Install from PyPi @@ -240,14 +240,18 @@ A.G. Cadavid, I. Montalban, A. Dalal, E. Solano, N.N. Hegade Efficient DCQO Algorithm within the Impulse Regime for Portfolio Optimization [arxiv:2308.15475](https://arxiv.org/abs/2308.15475) +A. Anand, K. Brown +Hamiltonians, groups, graphs and ansätze +[arxiv:2312.17146](https://arxiv.org/abs/2312.17146) + +P.W.K. Jensen, E.R. Kjellgren, P. Reinholdt, K.M. Ziems, S. Coriani, J. Kongsted, S. Sauer +Quantum Equation of Motion with Orbital Optimization for Computing Molecular Properties in Near-Term Quantum Computing +[arxiv:2312.12386](https://arxiv.org/abs/2312.12386) + Let us know, if you want your research project and/or tutorial to be included in this list! # Dependencies -Support for additional optimizers or quantum backends can be activated by intalling them in your environment. -Tequila will then detect them automatically. -Currently those are: [Phoenics](https://github.com/aspuru-guzik-group/phoenics) - and [GPyOpt](https://sheffieldml.github.io/GPyOpt/). -Quantum backends are treated in the same way. +Support for specific backends (quantum simulators, optimizers, quantum chemistry) can be activated by intalling them in your environment. ## Quantum Backends Currently supported @@ -267,7 +271,7 @@ Currently supported ### [Psi4](https://github.com/psi4/psi4). In a conda environment this can be installed with ```bash -conda install psi4 -c psi4 +conda install psi4 -c conda-forge ``` Here is a small [tutorial](https://nbviewer.org/github/tequilahub/tequila-tutorials/blob/main/chemistry/ChemistryModule.ipynb) that illustrates the usage. @@ -381,5 +385,3 @@ Tequila runs on Mac OSX. You might get in trouble with installing qulacs since it currently does not work with Apple's clang compiler. You need to install latest GNU compile (at least gcc-7 and g++7) and set them as default before installing qulacs over pip. -## Qibo and GPyOpt -Currently you can't use Qibo and GPyOpt within the same environment. diff --git a/requirements_gpyopt.txt b/requirements_gpyopt.txt deleted file mode 100644 index ac7fc7b0..00000000 --- a/requirements_gpyopt.txt +++ /dev/null @@ -1,17 +0,0 @@ -### requirements for gpyopt optimizer -cycler>=0.10.0 -decorator>=4.0.10 -numpy>=1.11.2 -six>=1.10.0 -python-dateutil>=2.6.0 -paramz>=0.7.0 -GPy>=1.8 -matplotlib>=1.5.3 -pyparsing>=2.1.10 -pytz>=2016.7 -scipy>=0.18.1 -mock>=2.0.0 -PyDOE >= 0.3.0 -sobol_seq >=0.1 -emcee>=2.2.1 -gpyopt diff --git a/src/tequila/circuit/gates.py b/src/tequila/circuit/gates.py index 3662e607..f39a9cb0 100644 --- a/src/tequila/circuit/gates.py +++ b/src/tequila/circuit/gates.py @@ -940,7 +940,7 @@ def QubitExcitation(angle: typing.Union[numbers.Real, Variable, typing.Hashable] except: raise Exception("QubitExcitation: Needs an even number of targets") - return QCircuit.wrap_gate(QubitExcitationImpl(angle=angle, target=target, assume_real=assume_real, compile_options=compile_options)) + return QCircuit.wrap_gate(QubitExcitationImpl(angle=angle, target=target, assume_real=assume_real, compile_options=compile_options, control=control)) """ diff --git a/src/tequila/circuit/qasm.py b/src/tequila/circuit/qasm.py index 65878f76..16a35b57 100644 --- a/src/tequila/circuit/qasm.py +++ b/src/tequila/circuit/qasm.py @@ -313,19 +313,33 @@ def parse_command(command: str, custom_gates_map: Dict[str, QCircuit], qregister return apply_custom_gate(custom_circuit=custom_circuit, qregisters_values=qregisters_values) if name in ("x", "y", "z", "h", "cx", "cy", "cz", "ch"): - return QCircuit.wrap_gate(gates.impl.QGateImpl(name=(name[1] if name[0] == 'c' else name).upper(), - control=get_qregister(args[0], qregisters) if name[0] == 'c' else None, - target=get_qregister(args[1 if name[0] == 'c' else 0], qregisters))) + target = get_qregister(args[0], qregisters) + control = None + if name[0].lower() == 'c': + control = get_qregister(args[0], qregisters) + target = get_qregister(args[1], qregisters) + name = name[1] + G = getattr(gates, name.upper()) + return G(control=control, target=target) + if name in ("ccx", "ccy", "ccz"): - return QCircuit.wrap_gate(gates.impl.QGateImpl(name=name.upper()[2], - control=[get_qregister(args[0], qregisters), get_qregister(args[1], qregisters)], - target=get_qregister(args[2], qregisters))) + G = getattr(gates, name[2].upper()) + control = [get_qregister(args[0], qregisters), get_qregister(args[1], qregisters)] + target = get_qregister(args[2], qregisters) + return G(control=control, target=target) + if name.startswith("rx(") or name.startswith("ry(") or name.startswith("rz(") or \ name.startswith("crx(") or name.startswith("cry(") or name.startswith("crz("): - return QCircuit.wrap_gate(gates.impl.RotationGateImpl(axis=name[2 if name[0] == 'c' else 1], - angle=get_angle(name)[0], - control=get_qregister(args[0], qregisters) if name[0] == 'c' else None, - target=get_qregister(args[1 if name[0] == 'c' else 0], qregisters))) + angle = get_angle(name)[0] + i = name.find('(') + name = name[0:i] + name = name.upper() + name = [x for x in name] + name[-1] = name[-1].lower() + name = "".join(name) + G = getattr(gates, name) + return G(angle=angle,control=get_qregister(args[0], qregisters) if name[0] == 'C' else None,target=get_qregister(args[1 if name[0] == 'C' else 0], qregisters)) + if name.startswith("U("): angles = get_angle(name) return gates.U(theta=angles[0], phi=angles[1], lambd=angles[2], @@ -362,7 +376,7 @@ def parse_command(command: str, custom_gates_map: Dict[str, QCircuit], qregister control=get_qregister(args[0], qregisters), target=get_qregister(args[1], qregisters)) if name in ("s", "t", "sdg", "tdg"): - g = gates.Phase(pi / (2 if name.startswith("s") else 4), + g = gates.Phase(angle=pi / (2 if name.startswith("s") else 4), control=None, target=get_qregister(args[0], qregisters)) if name.find("dg") != -1: diff --git a/src/tequila/circuit/qpic.py b/src/tequila/circuit/qpic.py index 18faaffa..be2e394e 100644 --- a/src/tequila/circuit/qpic.py +++ b/src/tequila/circuit/qpic.py @@ -224,7 +224,7 @@ def export_to(circuit, 'always_use_generators': True, 'group_together': "BARRIER" } - elif not hasattr("style", "items"): + elif not hasattr(style, "items"): raise Exception( "style needs to be `tequila`, or `standard` or `generators` or a dictionary, you gave: {}".format( str(style))) diff --git a/src/tequila/grouping/binary_rep.py b/src/tequila/grouping/binary_rep.py index 6aea305d..ac15efd0 100644 --- a/src/tequila/grouping/binary_rep.py +++ b/src/tequila/grouping/binary_rep.py @@ -26,7 +26,7 @@ def init_from_qubit_hamiltonian(cls, hamiltonian: QubitHamiltonian, n_qubits=Non del Hof.terms[()] hamiltonian = QubitHamiltonian.from_openfermion(Hof) if n_qubits is None: - n_qubits = hamiltonian.n_qubits + n_qubits = max(hamiltonian.qubits)+1 binary_terms = [ BinaryPauliString( p.binary(n_qubits).binary, diff --git a/src/tequila/quantumchemistry/__init__.py b/src/tequila/quantumchemistry/__init__.py index f566522f..e04398c7 100644 --- a/src/tequila/quantumchemistry/__init__.py +++ b/src/tequila/quantumchemistry/__init__.py @@ -95,6 +95,8 @@ def Molecule(geometry: str = None, if backend is None: if basis_set is None or basis_set.lower() in ["madness", "mra", "pno"]: backend = "madness" + basis_set = "mra" + parameters.basis_set = basis_set if orbital_type is not None and orbital_type.lower() not in ["pno", "mra-pno"]: warnings.warn("only PNOs supported as orbital_type without basis set. Setting to pno - You gave={}".format(orbital_type), TequilaWarning) orbital_type = "pno" diff --git a/src/tequila/quantumchemistry/chemistry_tools.py b/src/tequila/quantumchemistry/chemistry_tools.py index ff1f3ff1..faa36843 100644 --- a/src/tequila/quantumchemistry/chemistry_tools.py +++ b/src/tequila/quantumchemistry/chemistry_tools.py @@ -2,12 +2,12 @@ import typing import warnings from dataclasses import dataclass - +from copy import deepcopy +from numbers import Real import numpy -from tequila import BitString, QCircuit, TequilaException +from tequila import BitString, QCircuit, TequilaException,Variable,compile_circuit from tequila.circuit import gates - try: from openfermion.ops.representations import get_active_space_integrals # needs openfermion 1.3 except ImportError as E: @@ -50,16 +50,132 @@ def __init__(self, generator, p0, transformation, indices=None, *args, **kwargs) self._name = "FermionicExcitation" self.transformation = transformation self.indices = indices - + if not hasattr(indices[0],"__len__"): + self.indices = [(indices[2 * i], indices[2 * i+1]) for i in range(len(indices) // 2)] + self.sign = self.format_excitation_variables(self.indices) + self.indices = self.format_excitation_indices(self.indices) def compile(self, *args, **kwargs): if self.is_convertable_to_qubit_excitation(): target = [] for x in self.indices: for y in x: target.append(y) - return gates.QubitExcitation(target=target, angle=-self.parameter, control=self.control) + return gates.QubitExcitation(target=target, angle=self.parameter, control=self.control) else: - return gates.Trotterized(generator=self.generator, control=self.control, angle=self.parameter, steps=1) + if self.transformation.lower().strip("_") == "jordanwigner": + return self.fermionic_excitation(angle=self.sign*self.parameter, indices=self.indices, control=self.control,opt=False) + else: + return gates.Trotterized(generator=self.generator, control=self.control, angle=self.parameter, steps=1) + def format_excitation_indices(self, idx): + """ + Consistent formatting of excitation indices + idx = [(p0,q0),(p1,q1),...,(pn,qn)] + sorted as: p0pair[0]: + sig *= -1 + for pair in range(len(idx)-1): + if idx[pair+1][0]>idx[pair][0]: + sig *= -1 + return sig + def cCRy(self, target: int, dcontrol: typing.Union[list, int], control: typing.Union[list, int], + angle: typing.Union[Real, Variable, typing.Hashable], case: int = 1) -> QCircuit: + ''' + Compilation of CRy as on https://doi.org/10.1103/PhysRevA.102.062612 + If not control passed, Ry returned + Parameters + ---------- + case: if 1 employs eq. 12 from the paper, if 0 eq. 13 + ''' + if control is not None and not len(control): + control = None + if isinstance(dcontrol, int): + dcontrol = [dcontrol] + if not len(dcontrol): + return compile_circuit(gates.Ry(angle=angle, target=target, control=control)) + else: + if isinstance(angle, str): + angle = Variable(angle) + U = QCircuit() + aux = dcontrol[0] + ctr = deepcopy(dcontrol) + ctr.pop(0) + if case: + U += self.cCRy(target=target, dcontrol=ctr, angle=angle / 2, case=1, control=control) + gates.H( + aux) + gates.CNOT(target, aux) + U += self.cCRy(target=target, dcontrol=ctr, angle=-angle / 2, case=0, control=control) + gates.CNOT( + target, aux) + gates.H(aux) + else: + U += gates.H(aux) + gates.CNOT(target, aux) + self.cCRy(target=target, dcontrol=ctr, angle=-angle / 2, + case=0, control=control) + U += gates.CNOT(target, aux) + gates.H(aux) + self.cCRy(target=target, dcontrol=ctr, angle=angle / 2, + case=1, control=control) + return U + + def fermionic_excitation(self, angle: typing.Union[Real, Variable, typing.Hashable], indices: typing.List, + control: typing.Union[int, typing.List] = None, opt: bool = True) -> QCircuit: + ''' + Excitation [(i,j),(k,l)],... compiled following https://doi.org/10.1103/PhysRevA.102.062612 + opt: whether to optimized CNOT H CNOT --> Rz Rz CNOT Rz + ''' + lto = [] + lfrom = [] + if isinstance(indices,tuple) and not hasattr(indices[0],"__len__"): + indices = [(indices[2 * i], indices[2 * i + 1]) for i in range(len(indices) // 2)] + for pair in indices: + lfrom.append(pair[0]) + lto.append(pair[1]) + Upair = QCircuit() + if isinstance(angle, str) or isinstance(angle, tuple): + angle = Variable(angle) + for i in range(len(lfrom) - 1): + Upair += gates.CNOT(lfrom[i + 1], lfrom[i]) + Upair += gates.CNOT(lto[i + 1], lto[i]) + Upair += gates.X(lto[i]) + gates.X(lfrom[i]) + Upair += gates.CNOT(lto[-1], lfrom[-1]) + crt = lfrom[::-1] + lto + Uladder = QCircuit() + pairs = lfrom + lto + pairs.sort() + orbs = [] + for o in range(len(pairs) // 2): + orbs += [*range(pairs[2 * o] + 1, pairs[2 * o + 1])] + if len(orbs): + for o in range(len(orbs) - 1): + Uladder += gates.CNOT(orbs[o], orbs[o + 1]) + Uladder += compile_circuit(gates.CZ(orbs[-1], lto[-1])) + crt.pop(-1) + if control is not None and (isinstance(control, int) or len(control) == 1): + if isinstance(control, int): + crt.append(control) + else: + crt = crt + control + control = [] + Ur = self.cCRy(target=lto[-1], dcontrol=crt, angle=angle, control=control) + Upair2 = Upair.dagger() + if opt: + Ur.gates.pop(-1) + Ur.gates.pop(-1) + Upair2.gates.pop(0) + Ur += gates.Rz(numpy.pi / 2, target=lto[-1]) + gates.Rz(-numpy.pi / 2, target=lfrom[-1]) + Ur += gates.CNOT(lto[-1], lfrom[-1]) + gates.Rz(numpy.pi / 2, target=lfrom[-1]) + gates.H(lfrom[-1]) + return Upair + Uladder + Ur + Uladder.dagger() + Upair2 def __str(self): if self.indices is not None: @@ -804,20 +920,18 @@ class IntegralManager: _one_body_integrals: numpy.ndarray = None _two_body_integrals: NBodyTensor = None _constant_term: float = None - _basis_type: str = "unknown" _basis_name: str = "unknown" _orbital_type: str = "unknown" # e.g. "HF", "PNO", "native" _orbital_coefficients: numpy.ndarray = None _active_space: ActiveSpaceData = None _orbitals: typing.List[OrbitalData] = None - def __init__(self, one_body_integrals, two_body_integrals, basis_type="custom", + def __init__(self, one_body_integrals, two_body_integrals, basis_name="unknown", orbital_type="unknown", constant_term=0.0, orbital_coefficients=None, active_space=None, overlap_integrals=None, orbitals=None, *args, **kwargs): self._one_body_integrals = one_body_integrals self._two_body_integrals = two_body_integrals self._constant_term = constant_term - self._basis_type = basis_type self._basis_name = basis_name self._orbital_type = orbital_type @@ -956,9 +1070,16 @@ def transform_to_native_orbitals(self): """ c = self.get_orthonormalized_orbital_coefficients() self.orbital_coefficients=c - self._orbital_type="orthonormalized-{}-basis".format(self._orbital_type) + self._orbital_type="orthonormalized-{}-basis".format(self._basis_name) + + def is_unitary(self, U): + if len(U.shape) != 2: return False + if U.shape[0] != U.shape[1]: return False + test = (U.conj().T).dot(U) - numpy.eye(U.shape[0]) + if not numpy.isclose(numpy.linalg.norm(test), 0.0): return False + return True - def transform_orbitals(self, U): + def transform_orbitals(self, U, name=None): """ Transform orbitals Parameters @@ -969,10 +1090,12 @@ def transform_orbitals(self, U): ------- updates the structure with new orbitals: c = cU """ - c = self.orbital_coefficients - c = numpy.einsum("ix, xj -> ij", c, U, optimize="greedy") - self.orbital_coefficients = c - self._orbital_type += "-transformed" + assert self.is_unitary(U) + self.orbital_coefficients = numpy.einsum("ix, xj -> ij", self.orbital_coefficients, U, optimize="greedy") + if name is None: + self._orbital_type += "-transformed" + else: + self._orbital_type = name def get_integrals(self, orbital_coefficients=None, ordering="openfermion", ignore_active_space=False, *args, **kwargs): """ @@ -1001,7 +1124,9 @@ def get_integrals(self, orbital_coefficients=None, ordering="openfermion", ignor active_integrals = get_active_space_integrals(one_body_integrals=h, two_body_integrals=g, occupied_indices=self._active_space.frozen_reference_orbitals, active_indices=self._active_space.active_orbitals) + c = active_integrals[0] + c + h = active_integrals[1] g = NBodyTensor(elems=active_integrals[2], ordering="openfermion") g.reorder(to=ordering) @@ -1069,14 +1194,16 @@ def __str__(self): result += str(x) + "\n" return result - def print_basis_info(self, *args, **kwargs) -> None: - print("{:15} : {}".format("basis_type", self._basis_type), *args, **kwargs) + def print_basis_info(self, print_coefficients=True, *args, **kwargs) -> None: print("{:15} : {}".format("basis_name", self._basis_name), *args, **kwargs) print("{:15} : {}".format("orbital_type", self._orbital_type), *args, **kwargs) - print("{:15} : {}".format("orthogonal", self.basis_is_orthogonal()), *args, **kwargs) - print("{:15} : {}".format("functions", self.one_body_integrals.shape[0]), *args, **kwargs) + print("{:15} : {}".format("orthogonal basis", self.basis_is_orthogonal()), *args, **kwargs) + print("{:15} : {}".format("basis functions", self.one_body_integrals.shape[0]), *args, **kwargs) + print("{:15} : {}".format("active orbitals", [o.idx_total for o in self.active_orbitals]), *args, **kwargs) print("{:15} : {}".format("reference", [x.idx_total for x in self.reference_orbitals]), *args, **kwargs) + if not print_coefficients: return + print("Current Orbitals", *args, **kwargs) for i,x in enumerate(self.orbitals): print(x, *args, **kwargs) diff --git a/src/tequila/quantumchemistry/encodings.py b/src/tequila/quantumchemistry/encodings.py index 12d5ed46..d31ae761 100644 --- a/src/tequila/quantumchemistry/encodings.py +++ b/src/tequila/quantumchemistry/encodings.py @@ -2,16 +2,21 @@ Collections of Fermion-to-Qubit encodings known to tequila Most are Interfaces to OpenFermion """ +import abc + +from tequila import TequilaException from tequila.circuit.circuit import QCircuit -from tequila.circuit.gates import X +from tequila.circuit.gates import X, CNOT from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian import openfermion +import numpy + def known_encodings(): # convenience for testing and I/O - encodings= { - "JordanWigner":JordanWigner, - "BravyiKitaev":BravyiKitaev, + encodings = { + "JordanWigner": JordanWigner, + "BravyiKitaev": BravyiKitaev, "BravyiKitaevFast": BravyiKitaevFast, "BravyiKitaevTree": BravyiKitaevTree, "TaperedBravyiKitaev": TaperedBravyKitaev @@ -22,14 +27,14 @@ def known_encodings(): "ReorderedBravyiKitaev": lambda **kwargs: BravyiKitaev(up_then_down=True, **kwargs), "ReorderedBravyiKitaevTree": lambda **kwargs: BravyiKitaevTree(up_then_down=True, **kwargs), } - return {k.replace("_","").replace("-","").upper():v for k,v in encodings.items()} + return {k.replace("_", "").replace("-", "").upper(): v for k, v in encodings.items()} -class EncodingBase: +class EncodingBase(metaclass=abc.ABCMeta): # true if the encoding is fully integrated # false: can only do special things (like building the Hamiltionian) # but is not consistent with UCC gate generation - _ucc_support=False + _ucc_support = False @property def supports_ucc(self): @@ -37,20 +42,20 @@ def supports_ucc(self): @property def name(self): - prefix="" + prefix = "" if self.up_then_down: - prefix="Reordered" + prefix = "Reordered" if hasattr(self, "_name"): - return prefix+self._name + return prefix + self._name else: - return prefix+type(self).__name__ + return prefix + type(self).__name__ def __init__(self, n_electrons, n_orbitals, up_then_down=False, *args, **kwargs): self.n_electrons = n_electrons self.n_orbitals = n_orbitals self.up_then_down = up_then_down - def __call__(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> QubitHamiltonian: + def __call__(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> QubitHamiltonian: """ :param fermion_operator: an openfermion FermionOperator @@ -58,7 +63,8 @@ def __call__(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs The openfermion QubitOperator of this class ecoding """ if self.up_then_down: - op = openfermion.reorder(operator=fermion_operator, order_function=openfermion.up_then_down, num_modes=2*self.n_orbitals) + op = openfermion.reorder(operator=fermion_operator, order_function=openfermion.up_then_down, + num_modes=2 * self.n_orbitals) else: op = fermion_operator @@ -73,18 +79,18 @@ def up(self, i): if self.up_then_down: return i else: - return 2*i + return 2 * i def down(self, i): if self.up_then_down: - return i+self.n_orbitals + return i + self.n_orbitals else: - return 2*i+1 + return 2 * i + 1 - def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: + def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: raise Exception("{}::do_transform: called base class".format(type(self).__name__)) - def map_state(self, state:list, *args, **kwargs) -> list: + def map_state(self, state: list, *args, **kwargs) -> list: """ Expects a state in spin-orbital ordering Returns the corresponding qubit state in the class encoding @@ -112,7 +118,7 @@ def map_state(self, state:list, *args, **kwargs) -> list: # default is a lazy workaround, but it workds n_qubits = 2 * self.n_orbitals - spin_orbitals = sorted([i for i,x in enumerate(state) if int(x)==1]) + spin_orbitals = sorted([i for i, x in enumerate(state) if int(x) == 1]) string = "1.0 [" for i in spin_orbitals: @@ -128,26 +134,82 @@ def map_state(self, state:list, *args, **kwargs) -> list: key = list(wfn.keys())[0].array return key - def hcb_to_me(self, *args, **kwargs): - return None + @abc.abstractmethod + def me_to_jw(self) -> QCircuit: + """ + This method needs to be implemented to enable default conversions via Jordan-Wigner + """ + pass + + # independent conversion methods, these are used for default conversions + # arXiv:1808.10402 IV. B. 2, Eq. 57 + # original: https://doi.org/10.1063/1.4768229 + def _jw_to_bk(self) -> QCircuit: + U = QCircuit() # Constructs empty circuit + + flipper = False + for i in range(self.n_orbitals * 2): + # even qubits only hold their own value + if i % 2 == 0: + continue + + # sum always includes the last qubit + U += CNOT(control=i - 1, target=i) + + # every second odd qubit ties together with the last odd qubit + if flipper: + U += CNOT(control=i - 2, target=i) + + flipper = not flipper + + # we have now created the 4x4 blocks on the diagonal of this operators matrix + + # every power of 2 connects to the last power of 2 + # this corresponds to the last row in the recursive definitions being all 1s + x = numpy.log2(i + 1) + if x.is_integer() and x >= 3: + x = int(x) + U += CNOT(control=2 ** (x - 1) - 1, target=i) + + return U + + def _hcb_to_jw(self): + U = QCircuit() + for i in range(self.n_orbitals): + U += X(target=self.down(i), control=self.up(i)) + return U + + # Convenience Methods + def jw_to_me(self) -> QCircuit: + return self.me_to_jw().dagger() + + def me_to_bk(self) -> QCircuit: + return self.me_to_jw() + self._jw_to_bk() + + def bk_to_me(self) -> QCircuit: + return self.me_to_bk().dagger() + + def hcb_to_me(self) -> QCircuit: + return self._hcb_to_jw() + self.jw_to_me() def __str__(self): return type(self).__name__ + class JordanWigner(EncodingBase): """ OpenFermion::jordan_wigner """ - _ucc_support=True + _ucc_support = True - def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: + def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: return openfermion.jordan_wigner(fermion_operator, *args, **kwargs) - def map_state(self, state:list, *args, **kwargs): - state = state + [0]*(self.n_orbitals-len(state)) - result = [0]*len(state) + def map_state(self, state: list, *args, **kwargs): + state = state + [0] * (self.n_orbitals - len(state)) + result = [0] * len(state) if self.up_then_down: - return [state[2*i] for i in range(self.n_orbitals)] + [state[2*i+1] for i in range(self.n_orbitals)] + return [state[2 * i] for i in range(self.n_orbitals)] + [state[2 * i + 1] for i in range(self.n_orbitals)] else: return state @@ -157,15 +219,34 @@ def hcb_to_me(self, *args, **kwargs): U += X(target=self.down(i), control=self.up(i)) return U + def me_to_jw(self) -> QCircuit: + return QCircuit() + + def jw_to_me(self) -> QCircuit: + return QCircuit() + + class BravyiKitaev(EncodingBase): """ Uses OpenFermion::bravyi_kitaev """ - _ucc_support=True + _ucc_support = True + + def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: + return openfermion.bravyi_kitaev(fermion_operator, n_qubits=self.n_orbitals * 2) + + def me_to_jw(self) -> QCircuit: + return self._jw_to_bk().dagger() - def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: - return openfermion.bravyi_kitaev(fermion_operator, n_qubits=self.n_orbitals*2) + def jw_to_me(self) -> QCircuit: + return self._jw_to_bk() + + def bk_to_me(self) -> QCircuit: + return QCircuit() + + def me_to_bk(self) -> QCircuit: + return QCircuit() class BravyiKitaevTree(EncodingBase): @@ -173,28 +254,37 @@ class BravyiKitaevTree(EncodingBase): Uses OpenFermion::bravyi_kitaev_tree """ - _ucc_support=True + _ucc_support = True + + def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: + return openfermion.bravyi_kitaev_tree(fermion_operator, n_qubits=self.n_orbitals * 2) + + def me_to_jw(self) -> QCircuit: + raise TequilaException("{}::me_to_jw: unimplemented".format(type(self).__name__)) - def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: - return openfermion.bravyi_kitaev_tree(fermion_operator, n_qubits=self.n_orbitals*2) class BravyiKitaevFast(EncodingBase): """ Uses OpenFermion::bravyi_kitaev_tree """ - _ucc_support=False + _ucc_support = False - def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: + def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: n_qubits = openfermion.count_qubits(fermion_operator) - if n_qubits != self.n_orbitals*2: - raise Exception("BravyiKitaevFast transformation currently only possible for full Hamiltonians (no UCC generators).\nfermion_operator was {}".format(fermion_operator)) + if n_qubits != self.n_orbitals * 2: + raise Exception( + "BravyiKitaevFast transformation currently only possible for full Hamiltonians (no UCC generators).\nfermion_operator was {}".format( + fermion_operator)) op = openfermion.get_interaction_operator(fermion_operator) return openfermion.bravyi_kitaev_fast(op) -class TaperedBravyKitaev(EncodingBase): + def me_to_jw(self) -> QCircuit: + raise TequilaException("{}::me_to_jw: unimplemented".format(type(self).__name__)) - _ucc_support=False + +class TaperedBravyKitaev(EncodingBase): + _ucc_support = False """ Uses OpenFermion::symmetry_conserving_bravyi_kitaev (tapered bravyi_kitaev_tree arxiv:1701.07072) @@ -202,6 +292,7 @@ class TaperedBravyKitaev(EncodingBase): See OpenFermion Documentation for more Does not work for UCC generators yet """ + def __init__(self, n_electrons, n_orbitals, active_fermions=None, active_orbitals=None, *args, **kwargs): if active_fermions is None: self.active_fermions = n_electrons @@ -209,7 +300,7 @@ def __init__(self, n_electrons, n_orbitals, active_fermions=None, active_orbital self.active_fermions = active_fermions if active_orbitals is None: - self.active_orbitals = n_orbitals*2 # in openfermion those are spin-orbitals + self.active_orbitals = n_orbitals * 2 # in openfermion those are spin-orbitals else: self.active_orbitals = active_orbitals @@ -217,17 +308,20 @@ def __init__(self, n_electrons, n_orbitals, active_fermions=None, active_orbital raise Exception("Don't pass up_then_down argument to {}, it can't be changed".format(type(self).__name__)) super().__init__(n_orbitals=n_orbitals, n_electrons=n_electrons, up_then_down=False, *args, **kwargs) - def do_transform(self, fermion_operator:openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: - if openfermion.count_qubits(fermion_operator) != self.n_orbitals*2: + def do_transform(self, fermion_operator: openfermion.FermionOperator, *args, **kwargs) -> openfermion.QubitOperator: + if openfermion.count_qubits(fermion_operator) != self.n_orbitals * 2: raise Exception("TaperedBravyiKitaev not ready for UCC generators yet") - return openfermion.symmetry_conserving_bravyi_kitaev(fermion_operator, active_orbitals=self.active_orbitals, active_fermions=self.active_fermions) + return openfermion.symmetry_conserving_bravyi_kitaev(fermion_operator, active_orbitals=self.active_orbitals, + active_fermions=self.active_fermions) - def map_state(self, state:list, *args, **kwargs): - non_tapered_trafo = BravyiKitaevTree(up_then_down=True, n_electrons=self.n_electrons, n_orbitals=self.n_orbitals) + def map_state(self, state: list, *args, **kwargs): + non_tapered_trafo = BravyiKitaevTree(up_then_down=True, n_electrons=self.n_electrons, + n_orbitals=self.n_orbitals) key = non_tapered_trafo.map_state(state=state, *args, **kwargs) - n_qubits = self.n_orbitals*2 + n_qubits = self.n_orbitals * 2 active_qubits = [i for i in range(n_qubits) if i not in [n_qubits - 1, n_qubits // 2 - 1]] key = [key[i] for i in active_qubits] return key - + def me_to_jw(self) -> QCircuit: + raise TequilaException("{}::me_to_jw: unimplemented".format(type(self).__name__)) diff --git a/src/tequila/quantumchemistry/madness_interface.py b/src/tequila/quantumchemistry/madness_interface.py index 9273046c..7b291760 100644 --- a/src/tequila/quantumchemistry/madness_interface.py +++ b/src/tequila/quantumchemistry/madness_interface.py @@ -420,7 +420,10 @@ def make_upccgsd_ansatz(self, name="UpCCGSD", label=None, direct_compiling=None, """ # check if the used qubit encoding has a hcb transformation - have_hcb_trafo = self.transformation.hcb_to_me() is not None + try: + have_hcb_trafo = self.transformation.hcb_to_me() is not None + except: + have_hcb_trafo = False name = name.upper() # Default Method diff --git a/src/tequila/quantumchemistry/orbital_optimizer.py b/src/tequila/quantumchemistry/orbital_optimizer.py index a7bca8d7..5972134a 100644 --- a/src/tequila/quantumchemistry/orbital_optimizer.py +++ b/src/tequila/quantumchemistry/orbital_optimizer.py @@ -37,7 +37,7 @@ def __call__(self, local_data, *args, **kwargs): self.iterations += 1 def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=None, silent=False, - vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, molecule_factory=None, molecule_arguments=None, *args, **kwargs): + vqe_solver_arguments=None, initial_guess=None, return_mcscf=False, use_hcb=False, molecule_factory=None, molecule_arguments=None, restrict_to_active_space=True, *args, **kwargs): """ Parameters @@ -78,7 +78,12 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N if pyscf_arguments is None: pyscf_arguments = {"max_cycle_macro": 10, "max_cycle_micro": 3} no = molecule.n_orbitals - pyscf_molecule = QuantumChemistryPySCF.from_tequila(molecule=molecule, transformation=molecule.transformation) + + if not isinstance(molecule, QuantumChemistryPySCF): + pyscf_molecule = QuantumChemistryPySCF.from_tequila(molecule=molecule, transformation=molecule.transformation) + else: + pyscf_molecule = molecule + mf = pyscf_molecule._get_hf() result=OptimizeOrbitalsResult() mc = mcscf.CASSCF(mf, pyscf_molecule.n_orbitals, pyscf_molecule.n_electrons) @@ -140,10 +145,11 @@ def optimize_orbitals(molecule, circuit=None, vqe_solver=None, pyscf_arguments=N mc.kernel() # make new molecule - transformed_molecule = pyscf_molecule.transform_orbitals(orbital_coefficients=mc.mo_coeff) + mo_coeff = mc.mo_coeff + transformed_molecule = pyscf_molecule.transform_orbitals(orbital_coefficients=mo_coeff, name="optimized") result.molecule=transformed_molecule result.old_molecule=molecule - result.mo_coeff=mc.mo_coeff + result.mo_coeff=mo_coeff result.energy=mc.e_tot if return_mcscf: diff --git a/src/tequila/quantumchemistry/pyscf_interface.py b/src/tequila/quantumchemistry/pyscf_interface.py index bed3cc5d..69b681f9 100644 --- a/src/tequila/quantumchemistry/pyscf_interface.py +++ b/src/tequila/quantumchemistry/pyscf_interface.py @@ -75,6 +75,7 @@ def __init__(self, parameters: ParametersQC, kwargs["two_body_integrals"] = g_ao kwargs["one_body_integrals"] = h_ao kwargs["orbital_coefficients"] = mo_coeff + kwargs["orbital_type"] = "hf" if "nuclear_repulsion" not in kwargs: kwargs["nuclear_repulsion"] = mol.energy_nuc() diff --git a/src/tequila/quantumchemistry/qc_base.py b/src/tequila/quantumchemistry/qc_base.py index 9689b042..7364a58e 100644 --- a/src/tequila/quantumchemistry/qc_base.py +++ b/src/tequila/quantumchemistry/qc_base.py @@ -17,7 +17,7 @@ import typing, numpy, numbers from itertools import product - +import tequila.grouping.fermionic_functions as ff try: @@ -32,8 +32,7 @@ except Exception as E: raise Exception("{}\nIssue with Tequila Chemistry: Please update openfermion".format(str(E))) import warnings - - +OPTIMIZED_ORDERING = "Optimized" class QuantumChemistryBase: """ Base Class for tequila chemistry functionality @@ -94,7 +93,7 @@ def __init__(self, parameters: ParametersQC, else: self.integral_manager = self.initialize_integral_manager(active_orbitals=active_orbitals, reference_orbitals=reference_orbitals, - orbitals=orbitals, frozen_orbitals=frozen_orbitals, orbital_type=orbital_type, *args, + orbitals=orbitals, frozen_orbitals=frozen_orbitals, orbital_type=orbital_type, basis_name=self.parameters.basis_set, *args, **kwargs) if orbital_type is not None and orbital_type.lower() == "native": @@ -109,15 +108,28 @@ def __init__(self, parameters: ParametersQC, @classmethod def from_tequila(cls, molecule, transformation=None, *args, **kwargs): - c, h1, h2 = molecule.get_integrals() + c = molecule.integral_manager.constant_term + h1 = molecule.integral_manager.one_body_integrals + h2 = molecule.integral_manager.two_body_integrals + S = molecule.integral_manager.overlap_integrals + if "active_orbitals" not in kwargs: + active_orbitals = [o.idx_total for o in molecule.integral_manager.active_orbitals] + else: + active_orbitals = kwargs["active_orbitals"] + kwargs.pop("active_orbitals") if transformation is None: transformation = molecule.transformation + parameters = molecule.parameters return cls(nuclear_repulsion=c, one_body_integrals=h1, two_body_integrals=h2, - n_electrons=molecule.n_electrons, + overlap_integrals = S, + orbital_coefficients = molecule.integral_manager.orbital_coefficients, + active_orbitals= active_orbitals, transformation=transformation, - parameters=molecule.parameters, *args, **kwargs) + orbital_type=molecule.integral_manager._orbital_type, + parameters=parameters, + reference_orbitals= molecule.integral_manager.active_space.reference_orbitals,*args, **kwargs) def supports_ucc(self): """ @@ -433,10 +445,14 @@ def make_excitation_gate(self, indices, angle, control=None, assume_real=True, * generator = self.make_excitation_generator(indices=indices, remove_constant_term=control is None) p0 = self.make_excitation_generator(indices=indices, form="P0", remove_constant_term=control is None) - + if self.transformation.up_then_down: + idx = [] + for pair in indices: + idx.append((pair[0]//2+(pair[0]%2)*self.n_orbitals,pair[1]//2+(pair[1]%2)*self.n_orbitals)) + else:idx = indices return QCircuit.wrap_gate( FermionicGateImpl(angle=angle, generator=generator, p0=p0, - transformation=type(self.transformation).__name__.lower(), indices=indices, + transformation=type(self.transformation).__name__.lower(), indices=idx, assume_real=assume_real, control=control, **kwargs)) @@ -543,11 +559,13 @@ def initialize_integral_manager(self, *args, **kwargs): return manager - def transform_orbitals(self, orbital_coefficients, *args, **kwargs): + def transform_orbitals(self, orbital_coefficients, ignore_active_space=False, name=None, *args, **kwargs): """ Parameters ---------- - orbital_coefficients: second index is new orbital indes, first is old orbital index (summed over) + orbital_coefficients: second index is new orbital indes, first is old orbital index (summed over), indices are assumed to be defined on the active space + ignore_active_space: if true orbital_coefficients are not assumed to be given in the active space + name: str, name the new orbitals args kwargs @@ -556,9 +574,20 @@ def transform_orbitals(self, orbital_coefficients, *args, **kwargs): New molecule with transformed orbitals """ + U = numpy.eye(self.integral_manager.orbital_coefficients.shape[0]) + # mo_coeff by default only acts on the active space + active_indices = [o.idx_total for o in self.integral_manager.active_orbitals] + + if ignore_active_space: + U = orbital_coefficients + else: + for kk,k in enumerate(active_indices): + for ll,l in enumerate(active_indices): + U[k][l] = orbital_coefficients[kk][ll] + # can not be an instance of a specific backend (otherwise we get inconsistencies with classical methods in the backend) integral_manager = copy.deepcopy(self.integral_manager) - integral_manager.transform_orbitals(U=orbital_coefficients) + integral_manager.transform_orbitals(U=U, name=name) result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, transformation=self.transformation) return result @@ -583,7 +612,7 @@ def use_native_orbitals(self, inplace=False): else: integral_manager = copy.deepcopy(self.integral_manager) integral_manager.transform_to_native_orbitals() - result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, orbital_type="native", transformation=self.transformation) + result = QuantumChemistryBase(parameters=self.parameters, integral_manager=integral_manager, transformation=self.transformation) return result @@ -645,6 +674,68 @@ def n_electrons(self) -> int: """ return 2 * len(self.integral_manager.active_reference_orbitals) + def make_annihilation_op(self, orbital, coefficient=1.0): + """ + Compute annihilation operator on spin-orbital in qubit representation + Spin-orbital order is always (up,down,up,down,...) + """ + assert orbital<=self.n_orbitals*2 + aop = openfermion.ops.FermionOperator(f'{orbital}', coefficient) + return self.transformation(aop) + + def make_creation_op(self, orbital, coefficient=1.0): + """ + Compute creation operator on spin-orbital in qubit representation + Spin-orbital order is always (up,down,up,down,...) + """ + assert orbital<=self.n_orbitals*2 + cop = openfermion.ops.FermionOperator(f'{orbital}^', coefficient) + return self.transformation(cop) + + def make_number_op(self, orbital): + """ + Compute number operator on spin-orbital in qubit representation + Spin-orbital order is always (up,down,up,down,...) + """ + num_op = self.make_creation_op(orbital) * self.make_annihilation_op(orbital) + return num_op + + def make_sz_op(self): + """ + Compute the spin_z operator of the molecule in qubit representation + """ + sz = QubitHamiltonian() + for i in range(0, self.n_orbitals * 2, 2): + one = 0.5 * self.make_creation_op(i) * self.make_annihilation_op(i) + two = 0.5 * self.make_creation_op(i+1) * self.make_annihilation_op(i+1) + sz += (one - two) + return sz + + def make_sp_op(self): + """ + Compute the spin+ operator of the molecule in qubit representation + """ + sp = QubitHamiltonian() + for i in range(self.n_orbitals): + sp += self.make_creation_op(i*2) * self.make_annihilation_op(i*2 + 1) + return sp + + def make_sm_op(self): + """ + Compute the spin- operator of the molecule in qubit representation + """ + sm = QubitHamiltonian() + for i in range(self.n_orbitals): + sm += self.make_creation_op(i*2 + 1) * self.make_annihilation_op(i*2) + return sm + + def make_s2_op(self): + """ + Compute the spin^2 operator of the molecule in qubit representation + """ + s2_op = self.make_sm_op() * self.make_sp_op() + self.make_sz_op() * (self.make_sz_op() + 1) + return s2_op + def make_hamiltonian(self, *args, **kwargs) -> QubitHamiltonian: """ Parameters @@ -805,13 +896,13 @@ def hcb_to_me(self, U=None, condensed=False): """ if U is None: U = QCircuit() - - # consistency - consistency = [x < self.n_orbitals for x in U.qubits] - if not all(consistency): - warnings.warn( - "hcb_to_me: given circuit is not defined on the first {} qubits. Is this a HCB circuit?".format( - self.n_orbitals)) + else: + ups = [self.transformation.up(i.idx) for i in self.orbitals] + consistency = [x in ups for x in U.qubits] + if not all(consistency): + warnings.warn( + "hcb_to_me: given circuit is not defined on all first {} qubits. Is this a HCB circuit?".format( + self.n_orbitals)) # map to alpha qubits if condensed: @@ -1156,7 +1247,13 @@ def make_upccgsd_ansatz(self, indices = self.make_upccgsd_indices(key=name) # check if the used qubit encoding has a hcb transformation - have_hcb_trafo = self.transformation.hcb_to_me() is not None + have_hcb_trafo = True + try: + if self.transformation.hcb_to_me() is None: + have_hcb_trafo = False + except: + have_hcb_trafo = False + # consistency checks for optimization if have_hcb_trafo and hcb_optimization is None and include_reference: @@ -1404,7 +1501,8 @@ def make_uccsd_ansatz(self, trotter_steps: int = 1, factor = 1.0 / trotter_steps for step in range(trotter_steps): for idx, angle in indices.items(): - UCCSD += self.make_excitation_gate(indices=idx, angle=factor * angle) + converted = [(idx[2 * i], idx[2 * i + 1]) for i in range(len(idx) // 2)] + UCCSD += self.make_excitation_gate(indices=converted, angle=factor * angle) if hasattr(initial_amplitudes, "lower") and initial_amplitudes.lower() == "mp2" and parametrized and add_singles: # mp2 has no singles, need to initialize them here (if not parametrized initializling as 0.0 makes no sense though) @@ -1638,7 +1736,8 @@ def rdm2(self): return None def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_free: bool = True, - get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False): + get_rdm1: bool = True, get_rdm2: bool = True, ordering="dirac", use_hcb: bool = False, + rdm_trafo: QubitHamiltonian = None, decompose=None): """ Computes the one- and two-particle reduced density matrices (rdm1 and rdm2) given a unitary U. This method uses the standard ordering in physics as denoted below. @@ -1666,6 +1765,9 @@ def compute_rdms(self, U: QCircuit = None, variables: Variables = None, spin_fre get_rdm1, get_rdm2 : Set whether either one or both rdm1, rdm2 should be computed. If both are needed at some point, it is recommended to compute them at once. + rdm_trafo : + The rdm operators can be transformed, e.g., a^dagger_i a_j -> U^dagger a^dagger_i a_j U, + where U represents the transformation. The default is set to None, implying that U equas the identity. Returns ------- @@ -1910,8 +2012,18 @@ def _build_2bdy_operators_hcb() -> list: # Transform operator lists to QubitHamiltonians if (not use_hcb): qops = [_get_qop_hermitian(op) for op in qops] + # Compute expected values - evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables) + if rdm_trafo is None: + if decompose is not None: + print("MANIPULATED") + X = decompose(H=qops, U=U) + evals = simulate(X, variables=variables) + else: + evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables) + else: + qops = [rdm_trafo.dagger()*qops[i]*rdm_trafo for i in range(len(qops))] + evals = simulate(ExpectationValue(H=qops, U=U, shape=[len(qops)]), variables=variables) # Assemble density matrices # If self._rdm1, self._rdm2 exist, reset them if they are of the other spin-type @@ -2044,6 +2156,58 @@ def perturbative_f12_correction(self, rdm1: numpy.ndarray = None, rdm2: numpy.nd n_ri=n_ri, external_info=external_info, **kwargs) return correction.compute() + def n_rotation(self, i, phi): + ''' + Creates a quantum circuit that applies a phase rotation based on phi to both components (up and down) of a given qubit. + + Parameters: + - i (int): The index of the qubit to which the rotation will be applied. + - phi (float): The rotation angle. The actual rotation applied will be multiplied with -2 for both components. + + Returns: + - QCircuit: A quantum circuit object containing the sequence of rotations applied to the up and down components of the specified qubit. + ''' + + # Generate number operators for the up and down components of the qubit. + n_up = self.make_number_op(2*i) + n_down = self.make_number_op(2*i+1) + + # Start a new circuit and apply rotations to each component. + circuit = gates.GeneralizedRotation(generator = n_up, angle=-2*phi) + circuit += gates.GeneralizedRotation(generator = n_down, angle=-2*phi) + return circuit + + def get_givens_circuit(self, unitary, tol = 1e-12, ordering = OPTIMIZED_ORDERING): + ''' + Constructs a quantum circuit from a given real unitary matrix using Givens rotations. + + This method decomposes a unitary matrix into a series of Givens and Rz (phase) rotations, + then constructs and returns a quantum circuit that implements this sequence of rotations. + + Parameters: + - unitary (numpy.array): A real unitary matrix representing the transformation to implement. + - tol (float): A tolerance threshold below which matrix elements are considered zero. + - ordering (list of tuples or 'Optimized'): Custom ordering of indices for Givens rotations or 'Optimized' to generate them automatically. + + Returns: + - QCircuit: A quantum circuit implementing the series of rotations decomposed from the unitary. + ''' + # Decompose the unitary matrix into Givens and phase (Rz) rotations. + theta_list, phi_list = get_givens_decomposition(unitary, tol, ordering) + + # Initialize an empty quantum circuit. + circuit = QCircuit() + + # Add all Rz (phase) rotations to the circuit. + for phi in phi_list: + circuit += self.n_rotation(phi[1], phi[0]) + + # Add all Givens rotations to the circuit. + for theta in reversed(theta_list): + circuit += self.UR(theta[1], theta[2], theta[0]*2) + + return circuit + def print_basis_info(self): return self.integral_manager.print_basis_info() @@ -2064,3 +2228,140 @@ def __str__(self) -> str: result += "\nmore information with: self.print_basis_info()\n" return result + +def givens_matrix(n, p, q, theta): + ''' + Construct a complex Givens rotation matrix of dimension n by theta between rows/columns p and q. + ''' + ''' + Generates a Givens rotation matrix of size n x n to rotate by angle theta in the (p, q) plane. This matrix can be complex + + Parameters: + - n (int): The size of the Givens rotation matrix. + - p (int): The first index for the rotation plane. + - q (int): The second index for the rotation plane. + - theta (float): The rotation angle. + + Returns: + - numpy.array: The Givens rotation matrix. + ''' + matrix = numpy.eye(n) # Matrix to hold complex numbers + cos_theta = numpy.cos(theta) + sin_theta = numpy.sin(theta) + + # Directly assign cosine and sine without complex phase adjustment + matrix[p, p] = cos_theta + matrix[q, q] = cos_theta + matrix[p, q] = sin_theta + matrix[q, p] = -sin_theta + + return matrix + +def get_givens_decomposition(unitary, tol = 1e-12, ordering = OPTIMIZED_ORDERING, return_diagonal = False): + ''' + Decomposes a real unitary matrix into Givens rotations (theta) and Rz rotations (phi). + + Parameters: + - unitary (numpy.array): A real unitary matrix to decompose. It cannot be complex. + - tol (float): Tolerance for considering matrix elements as zero. Elements with absolute value less than tol are treated as zero. + - ordering (list of tuples or 'Optimized'): Custom ordering of indices for Givens rotations or 'Optimized' to generate them automatically. + - return_diagonal (bool): If True, the function also returns the diagonal matrix as part of the output. + + Returns: + - list: A list of tuples, each representing a Givens rotation. Each tuple contains the rotation angle theta and indices (i,j) of the rotation. + - list: A list of tuples, each representing an Rz rotation. Each tuple contains the rotation angle phi and the index (i) of the rotation. + - numpy.array (optional): The diagonal matrix after applying all Givens rotations, returned if return_diagonal is True. + ''' + U = unitary # no need to copy as we don't modify the original + U[abs(U) < tol] = 0 # Zeroing out the small elements as per the tolerance level. + n = U.shape[0] + + # Determine optimized ordering if specified. + if ordering == OPTIMIZED_ORDERING: + ordering = ff.depth_eff_order_mf(n) + + theta_list = [] + phi_list = [] + + def calcTheta(U, c, r): + '''Calculate and apply the Givens rotation for a specific matrix element.''' + t = numpy.arctan2(-U[r,c], U[r-1,c]) + theta_list.append((t, r, r-1)) + g = givens_matrix(n,r,r-1,t) + U = numpy.dot(g, U) + + return U + + # Apply and store Givens rotations as per the given or computed ordering. + if ordering is None: + for c in range(n): + for r in range(n-1, c, -1): + U = calcTheta(U, c, r) + else: + for r, c in ordering: + U = calcTheta(U, c, r) + + # Calculating the Rz rotations based on the phases of the diagonal elements. + # For real elements this means a 180 degree shift, i.e. a sign change. + for i in range(n): + ph = numpy.angle(U[i,i]) + phi_list.append((ph, i)) + + # Filtering out rotations without significance. + theta_list_new = [] + for i, theta in enumerate(theta_list): + if abs(theta[0] % (2*numpy.pi)) > tol: + theta_list_new.append(theta) + + phi_list_new = [] + for i, phi in enumerate(phi_list): + if abs(phi[0]) > tol: + phi_list_new.append(phi) + + if return_diagonal: + # Optionally return the resulting diagonal + return theta_list_new, phi_list_new, U + else: + return theta_list_new, phi_list_new + +def reconstruct_matrix_from_givens(n, theta_list, phi_list, to_real_if_possible = True, tol = 1e-12): + ''' + Reconstructs a matrix from given Givens rotations and Rz diagonal rotations. + This function is effectively an inverse of get_givens_decomposition, and therefore only works with data in the same format as its output. + + Parameters: + - n (int): The size of the unitary matrix to be reconstructed. + - theta_list (list of tuples): Each tuple contains (angle, i, j) representing a Givens rotation of `angle` radians, applied to rows/columns `i` and `j`. + - phi_list (list of tuples): Each tuple contains (angle, i), representing an Rz rotation by `angle` radians applied to the `i`th diagonal element. + - to_real_if_possible (bool): If True, converts the matrix to real if its imaginary part is effectively zero. + - tol (float): The tolerance whether to swap a complex rotation for a sign change. + + Returns: + - numpy.ndarray: The reconstructed complex or real matrix, depending on the `to_real_if_possible` flag and matrix composition. + ''' + # Start with an identity matrix + reconstructed = numpy.eye(n, dtype=complex) + + # Apply Rz rotations for diagonal elements + for phi in phi_list: + angle, i = phi + # Directly apply a sign flip if the rotation angle is π + if numpy.isclose(angle, numpy.pi, atol=tol): + reconstructed[i, i] *= -1 + else: + reconstructed[i, i] *= numpy.exp(1j * angle) + + # Apply Givens rotations in reverse order + for theta in reversed(theta_list): + angle, i, j = theta + g = givens_matrix(n, i, j, angle) + reconstructed = numpy.dot(g.conj().T, reconstructed) # Transpose of Givens matrix applied to the left + + # Convert matrix to real if its imaginary part is negligible unless disabled via to_real_if_possible + if to_real_if_possible: + # Directly apply a sign flip if the rotation angle is π + if numpy.all(reconstructed.imag == 0): + # Convert to real by taking the real part + reconstructed = reconstructed.real + + return reconstructed diff --git a/src/tequila/simulators/simulator_qiskit.py b/src/tequila/simulators/simulator_qiskit.py index 769baa14..607a42a3 100755 --- a/src/tequila/simulators/simulator_qiskit.py +++ b/src/tequila/simulators/simulator_qiskit.py @@ -138,7 +138,7 @@ class BackendCircuitQiskit(BackendCircuit): } numbering = BitNumbering.LSB - + def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise=None, device=None, *args, **kwargs): """ @@ -170,7 +170,7 @@ def __init__(self, abstract_circuit: QCircuit, variables, qubit_map=None, noise= 'Rx': (lambda c: c.rx, lambda c: c.mcrx), 'Ry': (lambda c: c.ry, lambda c: c.mcry), 'Rz': (lambda c: c.rz, lambda c: c.mcrz), - 'Phase': (lambda c: c.u1, lambda c: c.cu1), + 'Phase': (lambda c: c.p, lambda c: c.cp), 'SWAP': (lambda c: c.swap, lambda c: c.cswap), } diff --git a/src/tequila/tools/random_generators.py b/src/tequila/tools/random_generators.py index 07f3a7cd..434022fc 100644 --- a/src/tequila/tools/random_generators.py +++ b/src/tequila/tools/random_generators.py @@ -2,6 +2,7 @@ from tequila.circuit import gates from tequila.circuit.circuit import QCircuit from tequila.hamiltonian.qubit_hamiltonian import QubitHamiltonian +from scipy.stats import unitary_group, ortho_group def make_random_circuit(n_qubits: int, rotation_gates: list=['rx', 'ry', 'rz'], n_rotations: int=None, enable_controls: bool=None) -> QCircuit: @@ -75,3 +76,19 @@ def make_random_hamiltonian(n_qubits: int , paulis: list=['X','Y','Z'], n_ps: in H = QubitHamiltonian(ham) return H + +def generate_random_unitary(size, complex = False): + ''' + Generates a random unitary (or furthermore orthogonal if complex is False) matrix of a specified size. + + Parameters: + - size (int): The size of the unitary matrix to be generated. + - complex (bool, optional): Whether the unitary should be complex. + + Returns: + - numpy.ndarray: A randomly generated unitary matrix. + ''' + if complex: + return unitary_group.rvs(size) + else: + return ortho_group.rvs(size) \ No newline at end of file diff --git a/src/tequila/version.py b/src/tequila/version.py index 2875226c..862db8be 100644 --- a/src/tequila/version.py +++ b/src/tequila/version.py @@ -1,2 +1,2 @@ -__version__ = "1.9.5dev" +__version__ = "1.9.6dev" __author__ = "Tequila Developers " diff --git a/tests/test_chemistry.py b/tests/test_chemistry.py index 86b2dc33..d24c6bfa 100644 --- a/tests/test_chemistry.py +++ b/tests/test_chemistry.py @@ -8,6 +8,8 @@ from tequila.objective import ExpectationValue from tequila.quantumchemistry.encodings import known_encodings from tequila.simulators.simulator_api import simulate +import tequila.quantumchemistry.qc_base as qcb +import tequila.tools.random_generators as rg HAS_PYSCF = "pyscf" in qc.INSTALLED_QCHEMISTRY_BACKENDS HAS_PSI4 = "psi4" in qc.INSTALLED_QCHEMISTRY_BACKENDS @@ -171,8 +173,6 @@ def test_ucc_singles_psi4(): def do_test_ucc(qc_interface, parameters, result, trafo, backend="qulacs"): # check examples for comments psi4_interface = qc_interface(parameters=parameters, transformation=trafo) - - hqc = psi4_interface.make_hamiltonian() amplitudes = psi4_interface.compute_ccsd_amplitudes() U = psi4_interface.make_uccsd_ansatz(trotter_steps=1, initial_amplitudes=amplitudes, include_reference_ansatz=True) variables = amplitudes.make_parameter_dictionary() @@ -368,9 +368,9 @@ def test_hamiltonian_reduction(backend): @pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found") @pytest.mark.parametrize("assume_real", [True, False]) -@pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "tapered_bravyi_kitaev"]) +@pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "reordered_jordan_wigner"]) def test_fermionic_gates(assume_real, trafo): - mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.7\nLi 0.0 0.0 0.0", basis_set="sto-3g") + mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.7\nLi 0.0 0.0 0.0", basis_set="sto-3g",transformation=trafo) U1 = mol.prepare_reference() U2 = mol.prepare_reference() variable_count = {} @@ -587,8 +587,14 @@ def test_spa_consistency(geometry, name, optimize, transformation): mol = tq.Molecule(geometry=geometry, basis_set="sto-3g", transformation=transformation).use_native_orbitals() - if mol.transformation.hcb_to_me() is None: + # test compares HCB and non-HCB SPA implementations + # conversion needs to be implemented for comparisson + # if not, HCB is not used in circuit optimization anyways + try: + mol.transformation.hcb_to_me() + except: return + # doesn't need to make physical sense for the test edges = [] i = 0 @@ -725,3 +731,55 @@ def test_orbital_optimization_hcb(geometry): assert numpy.isclose(opt1.energy,opt2.energy,atol=1.e-5) assert time1 < time2 assert (numpy.isclose(opt1.mo_coeff,opt2.mo_coeff, atol=1.e-5)).all() + +@pytest.mark.parametrize("transformation", ["JordanWigner", "ReorderedJordanWigner", "BravyiKitaev", "BravyiKitaevTree"]) +@pytest.mark.parametrize("size", [2, 8]) +def test_givens_on_molecule(size, transformation): + # dummy one-electron integrals + h = numpy.ones(shape=[size,size]) + # dummy two-electron integrals + g = numpy.ones(shape=[size, size, size, size]) + + U = rg.generate_random_unitary(size) + + # transformed integrals + th = (U.T.dot(h)).dot(U) + tg = numpy.einsum("ijkx, xl -> ijkl", g, U, optimize='greedy') + tg = numpy.einsum("ijxl, xk -> ijkl", tg, U, optimize='greedy') + tg = numpy.einsum("ixkl, xj -> ijkl", tg, U, optimize='greedy') + tg = numpy.einsum("xjkl, xi -> ijkl", tg, U, optimize='greedy') + + # original molecule/H + mol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=h, two_body_integrals=g, basis_set="dummy", transformation=transformation) + H = mol.make_hamiltonian() + # transformed molecule/H + tmol = tq.Molecule(geometry="He 0.0 0.0 0.0", nuclear_repulsion=0.0, one_body_integrals=th, two_body_integrals=tg,basis_set="dummy", transformation=transformation) + tH = tmol.make_hamiltonian() + + # transformation in qubit space (this corresponds to the U above) + UR = mol.get_givens_circuit(U) # Works! + + # test circuit + circuit = rg.make_random_circuit(size) + + # create expectation values and see if they are the same + E1 = tq.ExpectationValue(U=circuit, H=tH) + E2 = tq.ExpectationValue(U=circuit + UR, H=H) + + result1 = tq.simulate(E1) + result2 = tq.simulate(E2) + + assert numpy.isclose(result1, result2) + +@pytest.mark.parametrize("size", [2, 8]) +def test_givens_decomposition(size): + # generate random unitary + unitary = rg.generate_random_unitary(size) + + # decompose givens + theta_list, phi_list = qcb.get_givens_decomposition(unitary) + + # reconstruct original unitary from givens + reconstructed_matrix = qcb.reconstruct_matrix_from_givens(unitary.shape[0], theta_list, phi_list) + + assert numpy.allclose(unitary, reconstructed_matrix) diff --git a/tests/test_chemistry_madness.py b/tests/test_chemistry_madness.py index ae22ea7e..b398e26c 100644 --- a/tests/test_chemistry_madness.py +++ b/tests/test_chemistry_madness.py @@ -135,7 +135,7 @@ def test_madness_upccgsd(trafo): E = tq.ExpectationValue(H=H, U=U) assert (len(E.extract_variables()) == 2) variables = result.variables - if "bravyi" in trafo.lower(): + if "bravyikitaevtree" in trafo.lower(): # signs of angles change in BK compared to JW-like HCB variables = {k: -v for k, v in variables.items()} energy = tq.simulate(E, variables)