1# Licensed under the Apache License, Version 2.0 (the "License"); 2# you may not use this file except in compliance with the License. 3# You may obtain a copy of the License at 4# 5# http://www.apache.org/licenses/LICENSE-2.0 6# 7# Unless required by applicable law or agreed to in writing, software 8# distributed under the License is distributed on an "AS IS" BASIS, 9# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 10# See the License for the specific language governing permissions and 11# limitations under the License. 12 13import abc 14import itertools 15from typing import (cast, Dict, Optional, Sequence, Tuple, TYPE_CHECKING, Union) 16 17import numpy as np 18import scipy.linalg as la 19import sympy 20import cirq 21 22import openfermion.ops as ops 23from openfermion.transforms import get_interaction_operator 24from openfermion.linalg import jordan_wigner_sparse 25from openfermion.utils import hermitian_conjugated 26 27if TYPE_CHECKING: 28 import openfermion 29 30 31def _arg(x): 32 if x == 0: 33 return 0 34 if cirq.is_parameterized(x): 35 return sympy.arg(x) 36 return np.angle(x) 37 38 39def _canonicalize_weight(w): 40 if w == 0: 41 return (0, 0) 42 if cirq.is_parameterized(w): 43 return (cirq.PeriodicValue(abs(w), 2 * sympy.pi), sympy.arg(w)) 44 period = 2 * np.pi 45 return (np.round((w.real % period) if (w == np.real(w)) else 46 (abs(w) % period) * w / abs(w), 8), 0) 47 48 49def state_swap_eigen_component(x: str, y: str, sign: int = 1, angle: float = 0): 50 """The +/- eigen-component of the operation that swaps states x and y. 51 52 For example, state_swap_eigen_component('01', '10', ±1) with angle θ returns 53 ┌ ┐ 54 │0, 0, 0, 0│ 55 │0, 0.5, ±0.5 e^{-iθ}, 0│ 56 │0, ±0.5 e^{iθ}, 0.5, 0│ 57 │0, 0, 0, 0│ 58 └ ┘ 59 60 Args: 61 x: The first state to swap, as a bitstring. 62 y: The second state to swap, as a bitstring. Must have high index than 63 x. 64 sign: The sign of the off-diagonal elements (indicated by +/-1). 65 angle: The phase of the complex off-diagonal elements. Defaults to 0. 66 67 Returns: The eigen-component. 68 69 Raises: 70 ValueError: 71 * x and y have different lengths 72 * x or y contains a character other than '0' and '1' 73 * x and y are the same 74 * sign is not -1 or 1 75 TypeError: x or y is not a string 76 """ 77 if not (isinstance(x, str) and isinstance(y, str)): 78 raise TypeError('not (isinstance(x, str) and isinstance(y, str))') 79 if len(x) != len(y): 80 raise ValueError('len(x) != len(y)') 81 if set(x).union(y).difference('01'): 82 raise ValueError('Arguments must be 0-1 strings.') 83 if x == y: 84 raise ValueError('x == y') 85 if sign not in (-1, 1): 86 raise ValueError('sign not in (-1, 1)') 87 88 dim = 2**len(x) 89 i, j = int(x, 2), int(y, 2) 90 91 component = np.zeros((dim, dim), dtype=np.complex128) 92 component[i, i] = component[j, j] = 0.5 93 component[j, i] = sign * 0.5 * 1j**(angle * 2 / np.pi) 94 component[i, j] = sign * 0.5 * 1j**(-angle * 2 / np.pi) 95 return component 96 97 98def fermionic_simulation_gates_from_interaction_operator( 99 operator: 'openfermion.InteractionOperator'): 100 r""" 101 Given $H = \sum_{I \subset [n]} H_I$, returns gates 102 $\left\{G_I\right\} = \left\{e^{i H_I\right\}$. 103 104 Each term $H_I$ is the sum of all terms in $H$ that involve exactly the 105 orbitals $I$. 106 107 Args: 108 operator: The interaction operator ($H$). 109 110 Returns: A dict from tuples of mode indices to gates. 111 """ 112 n_qubits = operator.n_qubits 113 114 gates: Dict[Tuple[int, ...], cirq.Gate] = {} 115 116 if operator.constant: 117 gates[()] = operator.constant 118 for p in range(n_qubits): 119 coeff = operator.one_body_tensor[p, p] 120 if coeff: 121 gates[(p,)] = cirq.Z**(coeff / np.pi) 122 for modes in itertools.combinations(range(n_qubits), 2): 123 gate: Optional[InteractionOperatorFermionicGate] = ( 124 QuadraticFermionicSimulationGate.from_interaction_operator( 125 operator=operator, modes=modes)) 126 if gate: 127 gates[modes] = gate 128 for modes in itertools.combinations(range(n_qubits), 3): 129 gate = CubicFermionicSimulationGate.from_interaction_operator( 130 operator=operator, modes=modes) 131 if gate: 132 gates[modes] = gate 133 for modes in itertools.combinations(range(n_qubits), 4): 134 gate = QuarticFermionicSimulationGate.from_interaction_operator( 135 operator=operator, modes=modes) 136 if gate: 137 gates[modes] = gate 138 return gates 139 140 141def sum_of_interaction_operator_gate_generators( 142 n_modes: int, 143 gates: Dict[Tuple[int, ...], Union[float, cirq.Gate]], 144) -> 'openfermion.InteractionOperator': 145 """The interaction operator that is the sum of the generators of the 146 specified fermionic simulation gates. 147 148 The gates are specified as a dictionary whose items are (indices, gate), 149 where 150 * indices is a tuple of ints specifying the modes on which the gate 151 acts; 152 * gate is one of type 153 - float, which is interpreted as a constant, regardless of the 154 indices, 155 - cirq.ZPowGate, which is interpreted as a "linear" fermionic 156 simulation gate, 157 - openfermion.InteractionOperatorFermionicGate. 158 159 Args: 160 n_modes: The number of modes. 161 gates: The gates. 162 163 Returns: 164 The interaction operator. 165 """ 166 # assumes gate indices in JW order 167 operator = ops.InteractionOperator.zero(n_modes) 168 169 for indices, gate in gates.items(): 170 if not indices: 171 operator.constant += gate 172 elif isinstance(gate, cirq.ZPowGate): 173 coeff = gate._exponent * np.pi 174 operator.constant += gate._exponent * gate._global_shift * np.pi 175 operator.one_body_tensor[indices * 2] += coeff 176 elif isinstance(gate, InteractionOperatorFermionicGate): 177 gate.interaction_operator_generator(operator=operator, 178 modes=indices) 179 else: 180 raise TypeError(f'Gate type {gate} not supported.') 181 182 return operator 183 184 185@cirq.value_equality(approximate=True) 186class ParityPreservingFermionicGate(cirq.Gate, metaclass=abc.ABCMeta): 187 r"""The Jordan-Wigner transform of $\exp(-i H)$ for a fermionic 188 Hamiltonian $H$. 189 190 Each subclass corresponds to a set of generators $\{G_i\}$ 191 corresponding to the family of Hamiltonians $\sum_i w_i G_i + 192 \text{h.c.}$, where the weights $w_i \in \mathbb C$ are specified by 193 the instance. 194 195 The Jordan-Wigner mapping maps the fermionic modes $(0, \ldots, n - 196 1)$ to qubits $(0, \ldots, n - 1)$, respectively. 197 198 Each generator $G_i$ must be a linear combination of fermionic 199 monomials consisting of an even number of creation/annihilation operators. 200 This is so that the Jordan-Wigner transform acts only on the gate's qubits, 201 even when the fermionic modes are offset as part of a larger Jordan-Wigner 202 string. 203 """ 204 205 def __init__( 206 self, 207 weights: Optional[Tuple[complex, ...]] = None, 208 absorb_exponent: bool = False, 209 exponent: cirq.TParamVal = 1.0, 210 global_shift: float = 0.0, 211 ) -> None: 212 """A fermionic interaction. 213 214 Args: 215 weights: The weights of the terms in the Hamiltonian. 216 absorb_exponent: Whether to absorb the given exponent into the 217 weights. If true, the exponent of the return gate is `1`. 218 Defaults to `False`. 219 """ 220 if weights is None: 221 weights = (1.,) * self.num_weights() 222 self.weights = weights 223 224 self._exponent = exponent 225 self._global_shift = global_shift 226 self._canonical_exponent_cached = None 227 228 if absorb_exponent: 229 self.absorb_exponent_into_weights() 230 231 @staticmethod 232 @abc.abstractmethod 233 def fermion_generator_components() -> Tuple['openfermion.FermionOperator']: 234 r"""The FermionOperators $(G_i)_i$ such that the gate's fermionic 235 generator is $\sum_i w_i G_i + \text{h.c.}$ where $(w_i)_i$ 236 are the gate's weights.""" 237 238 @abc.abstractmethod 239 def fswap(self, i: int): 240 """Update the weights of the gate to effect conjugation by an FSWAP on 241 the i-th and (i+1)th qubits.""" 242 243 @classmethod 244 def num_weights(cls) -> int: 245 """The number of parameters (weights) in the generator.""" 246 return len(cls.fermion_generator_components()) 247 248 @property 249 def qubit_generator_matrix(self) -> np.ndarray: 250 """The matrix G such that the gate's unitary is exp(-i t G) with 251 exponent t.""" 252 return jordan_wigner_sparse(self.fermion_generator, 253 self.num_qubits()).toarray() 254 255 @property 256 def fermion_generator(self) -> 'openfermion.FermionOperator': 257 """The FermionOperator G such that the gate's unitary is exp(-i t G) 258 with exponent t using the Jordan-Wigner transformation.""" 259 half_generator = sum(( 260 w * G 261 for w, G in zip(self.weights, self.fermion_generator_components())), 262 ops.FermionOperator()) 263 return half_generator + hermitian_conjugated(half_generator) 264 265 def _diagram_exponent(self, 266 args: cirq.CircuitDiagramInfoArgs, 267 *, 268 ignore_global_phase: bool = True): 269 if not isinstance(self._exponent, (int, float)): 270 return self._exponent 271 result = float(self._exponent) 272 if args.precision is not None: 273 result = np.around(result, args.precision) 274 return result 275 276 @classmethod 277 def wire_symbol(cls, use_unicode: bool): 278 """The symbol to use in circuit diagrams.""" 279 return cls.__name__ 280 281 def _resolve_parameters_(self, resolver, recursive: bool = True): 282 resolved_weights = cirq.resolve_parameters(self.weights, 283 resolver, 284 recursive=recursive) 285 resolved_exponent = cirq.resolve_parameters(self._exponent, 286 resolver, 287 recursive=recursive) 288 resolved_global_shift = cirq.resolve_parameters(self._global_shift, 289 resolver, 290 recursive=recursive) 291 return type(self)(resolved_weights, 292 exponent=resolved_exponent, 293 global_shift=resolved_global_shift) 294 295 def _value_equality_values_(self): 296 return tuple( 297 _canonicalize_weight(w * self.exponent) 298 for w in list(self.weights) + [self._global_shift]) 299 300 def _is_parameterized_(self) -> bool: 301 return any( 302 cirq.is_parameterized(v) 303 for V in self._value_equality_values_() 304 for v in V) 305 306 def absorb_exponent_into_weights(self): 307 period = (2 * sympy.pi) if self._is_parameterized_() else 2 * (np.pi) 308 new_weights = [] 309 for weight in self.weights: 310 if not weight: 311 new_weights.append(weight) 312 continue 313 old_abs = abs(weight) 314 new_abs = (old_abs * self._exponent) % period 315 new_weights.append(weight * new_abs / old_abs) 316 self.weights = tuple(new_weights) 317 self._global_shift *= self._exponent 318 self._exponent = 1 319 320 def permute(self, init_pos: Sequence[int]): 321 """An in-place version of permuted.""" 322 I = range(self.num_qubits()) 323 if sorted(init_pos) != list(I): 324 raise ValueError(f'{init_pos} is not a permutation of {I}.') 325 curr_pos = list(init_pos) 326 for i in I: 327 for j in I[i % 2:-1:2]: 328 if curr_pos[j] > curr_pos[j + 1]: 329 self.fswap(j) 330 curr_pos[j:j + 2] = reversed(curr_pos[j:j + 2]) 331 assert curr_pos == list(I) 332 333 def permuted(self, init_pos: Sequence[int]): 334 """Returns a gate with the Jordan-Wigner ordering changed. 335 336 If the Jordan-Wigner ordering of the original gate is given by 337 init_pos, then the returned gate has Jordan-Wigner ordering 338 (0, ..., n - 1), where n is the number of qubits on which the gate acts. 339 340 Args: 341 init_pos: A permutation of (0, ..., n - 1). 342 """ 343 gate = self.__copy__() 344 gate.permute(init_pos) 345 return gate 346 347 def __copy__(self): 348 return type(self)(self.weights, 349 exponent=self.exponent, 350 global_shift=self._global_shift) 351 352 def _circuit_diagram_info_(self, args: cirq.CircuitDiagramInfoArgs 353 ) -> cirq.CircuitDiagramInfo: 354 wire_symbols = [self.wire_symbol(args.use_unicode_characters) 355 ] * self.num_qubits() 356 wire_symbols[0] += f'{tuple(self.weights)}' 357 exponent = self._diagram_exponent(args) 358 return cirq.CircuitDiagramInfo(wire_symbols=wire_symbols, 359 exponent=exponent) 360 361 362class InteractionOperatorFermionicGate(ParityPreservingFermionicGate): 363 r"""The Jordan-Wigner transform of $\exp(-i H)$ for a fermionic 364 Hamiltonian $H$, where $H$ is an interaction operator. 365 366 See openfermion.ParityPreservingFermionicGate and 367 openfermion.InteractionOperator for more details. 368 """ 369 370 @classmethod 371 @abc.abstractmethod 372 def from_interaction_operator( 373 cls, 374 *, 375 operator: 'openfermion.InteractionOperator', 376 modes: Optional[Sequence[int]] = None, 377 ) -> Optional['ParityPreservingFermionicGate']: 378 """Constructs the gate corresponding to the specified term in the 379 Hamiltonian.""" 380 381 def interaction_operator_generator( 382 self, 383 *, 384 operator: Optional['openfermion.InteractionOperator'] = None, 385 modes: Optional[Sequence[int]] = None 386 ) -> 'openfermion.InteractionOperator': 387 """Constructs the Hamiltonian corresponding to the gate's generator.""" 388 if modes is None: 389 modes = tuple(range(self.num_qubits())) 390 if operator is None: 391 n_modes = max(modes) + 1 392 operator = ops.InteractionOperator.zero(n_modes) 393 else: 394 n_modes = operator.n_qubits 395 fermion_operator = self.fermion_generator 396 return get_interaction_operator(fermion_operator, n_qubits=n_modes) 397 398 399class QuadraticFermionicSimulationGate(InteractionOperatorFermionicGate, 400 cirq.InterchangeableQubitsGate, 401 cirq.EigenGate): 402 r"""``(w0 |10⟩⟨01| + h.c.) + w1 |11⟩⟨11|`` interaction. 403 404 With weights $(w_0, w_1)$ and exponent $t$, this gate's matrix 405 is defined as 406 407 $$ 408 e^{-i t H}, 409 $$ 410 411 where 412 413 $$ 414 H = \left(w_0 \left| 10 \right\rangle\left\langle 01 \right| + 415 \text{h.c.}\right) - 416 w_1 \left| 11 \right\rangle \left\langle 11 \right|. 417 $$ 418 419 This corresponds to the Jordan-Wigner transform of 420 421 $$ 422 H = (w_0 a^{\dagger}_i a_{i+1} + \text{h.c.}) + 423 w_1 a_{i}^{\dagger} a_{i+1}^{\dagger} a_{i} a_{i+1}, 424 $$ 425 426 where $a_i$ and $a_{i+1}$ are the annihilation operators for 427 the fermionic modes $i$ and $(i+1)$, respectively mapped to the 428 first and second qubits on which this gate acts. 429 430 Args: 431 weights: The weights of the terms in the Hamiltonian. 432 """ 433 434 @classmethod 435 def num_weights(cls): 436 return 2 437 438 def _num_qubits_(self) -> int: 439 return 2 440 441 def _decompose_(self, qubits): 442 r = 2 * abs(self.weights[0]) / np.pi 443 theta = _arg(self.weights[0]) / np.pi 444 yield cirq.Z(qubits[0])**-theta 445 yield cirq.ISwapPowGate(exponent=-r * self.exponent)(*qubits) 446 yield cirq.Z(qubits[0])**theta 447 yield cirq.CZPowGate(exponent=-self.weights[1] * self.exponent / 448 np.pi)(*qubits) 449 450 def _eigen_components(self): 451 components = [(0, np.diag([1, 0, 0, 0])), 452 (-self.weights[1] / np.pi, np.diag([0, 0, 0, 1]))] 453 r = abs(self.weights[0]) / np.pi 454 theta = 2 * _arg(self.weights[0]) / np.pi 455 for s in (-1, 1): 456 components.append( 457 (-s * r, 458 np.array([[0, 0, 0, 0], [0, 1, s * 1j**(-theta), 0], 459 [0, s * 1j**(theta), 1, 0], [0, 0, 0, 0]]) / 2)) 460 return components 461 462 def __repr__(self): 463 exponent_str = ('' if self.exponent == 1 else ', exponent=' + 464 cirq._compat.proper_repr(self.exponent)) 465 return ('openfermion.QuadraticFermionicSimulationGate(({}){})'.format( 466 ', '.join(cirq._compat.proper_repr(v) for v in self.weights), 467 exponent_str)) 468 469 @property 470 def qubit_generator_matrix(self): 471 generator = np.zeros((4, 4), dtype=np.complex128) 472 # w0 |10><01| + h.c. 473 generator[2, 1] = self.weights[0] 474 generator[1, 2] = self.weights[0].conjugate() 475 # w1 |11><11| 476 generator[3, 3] = self.weights[1] 477 return generator 478 479 @staticmethod 480 def fermion_generator_components(): 481 return ( 482 ops.FermionOperator(((0, 1), (1, 0))), 483 ops.FermionOperator(((0, 1), (0, 0), (1, 1), (1, 0)), 0.5), 484 ) 485 486 @classmethod 487 def wire_symbol(cls, use_unicode: bool): 488 return '↓↑' if use_unicode else 'a*a' 489 490 @classmethod 491 def from_interaction_operator( 492 cls, 493 *, 494 operator: 'openfermion.InteractionOperator', 495 modes: Optional[Sequence[int]] = None, 496 ) -> Optional['QuadraticFermionicSimulationGate']: 497 if modes is None: 498 modes = (0, 1) 499 500 p, q = modes 501 tunneling_coeff = operator.one_body_tensor[p, q] 502 interaction_coeff = (-operator.two_body_tensor[p, q, p, q] + 503 operator.two_body_tensor[q, p, p, q] + 504 operator.two_body_tensor[p, q, q, p] - 505 operator.two_body_tensor[q, p, q, p]) 506 weights: Tuple[complex, complex] = (-tunneling_coeff, 507 -interaction_coeff) 508 if any(weights): 509 return cls(weights) 510 return None 511 512 def interaction_operator_generator( 513 self, 514 *, 515 operator: Optional['openfermion.InteractionOperator'] = None, 516 modes: Optional[Sequence[int]] = None 517 ) -> 'openfermion.InteractionOperator': 518 if modes is None: 519 modes = (0, 1) 520 if operator is None: 521 n_modes = max(modes) + 1 522 operator = ops.InteractionOperator.zero(n_modes) 523 524 weights = tuple(w * self._exponent for w in self.weights) 525 operator.constant += self._exponent * self._global_shift 526 operator.one_body_tensor[modes] -= weights[0] 527 operator.one_body_tensor[modes[::-1]] -= weights[0].conjugate() 528 operator.two_body_tensor[tuple(modes) * 2] += weights[1] 529 return operator 530 531 def fswap(self, i: int = 0): 532 if i != 0: 533 raise ValueError(f'{i} != 0') 534 self.weights = (self.weights[0].conjugate(), self.weights[1]) 535 536 537class CubicFermionicSimulationGate(InteractionOperatorFermionicGate, 538 cirq.EigenGate): 539 r"""``w0|110⟩⟨101| + w1|110⟩⟨011| + w2|101⟩⟨011|`` + h.c. interaction. 540 541 With weights $(w_0, w_1, w_2)$ and exponent $t$, this gate's 542 matrix is defined as 543 544 $$ 545 e^{-i t H}, 546 $$ 547 548 where 549 550 $$ 551 H = \left(w_0 \left| 110 \right\rangle\left\langle 101 \right| + 552 \text{h.c.}\right) + 553 \left(w_1 \left| 110 \right\rangle\left\langle 011 \right| + 554 \text{h.c.}\right) + 555 \left(w_2 \left| 101 \right\rangle\left\langle 011 \right| + 556 \text{h.c.}\right) 557 $$ 558 559 This corresponds to the Jordan-Wigner transform of 560 561 $$ 562 H = -\left(w_0 a^{\dagger}_i a^{\dagger}_{i+1} a_{i} a_{i+2} + 563 \text{h.c.}\right) - 564 \left(w_1 a^{\dagger}_i a^{\dagger}_{i+1} a_{i+1} a_{i+2} + 565 \text{h.c.}\right) - 566 \left(w_2 a^{\dagger}_i a^{\dagger}_{i+2} a_{i+1} a_{i+2} + 567 \text{h.c.}\right), 568 $$ 569 570 where $a_i$, $a_{i+1}$, $a_{i+2}$ are the annihilation 571 operators for the fermionic modes $i$, $(i+1)$ $(i+2)$, 572 respectively mapped to the three qubits on which this gate acts. 573 574 Args: 575 weights: The weights of the terms in the Hamiltonian. 576 """ 577 578 @classmethod 579 def num_weights(cls): 580 return 3 581 582 def _num_qubits_(self) -> int: 583 return 3 584 585 @classmethod 586 def wire_symbol(cls, use_unicode: bool): 587 return '↕↓↑' if use_unicode else 'na*a' 588 589 def _eigen_components(self): 590 components = [(0, np.diag([1, 1, 1, 0, 1, 0, 0, 1]))] 591 nontrivial_part = np.zeros((3, 3), dtype=np.complex128) 592 for ij, w in zip([(1, 2), (0, 2), (0, 1)], self.weights): 593 nontrivial_part[ij] = w 594 nontrivial_part[ij[::-1]] = w.conjugate() 595 assert np.allclose(nontrivial_part, nontrivial_part.conj().T) 596 eig_vals, eig_vecs = np.linalg.eigh(nontrivial_part) 597 for eig_val, eig_vec in zip(eig_vals, eig_vecs.T): 598 exp_factor = -eig_val / np.pi 599 proj = np.zeros((8, 8), dtype=np.complex128) 600 nontrivial_indices = np.array([3, 5, 6], dtype=np.intp) 601 proj[nontrivial_indices[:, np.newaxis], nontrivial_indices] = ( 602 np.outer(eig_vec.conjugate(), eig_vec)) 603 components.append((exp_factor, proj)) 604 return components 605 606 def __repr__(self): 607 return ('openfermion.CubicFermionicSimulationGate(' + '({})'.format( 608 ' ,'.join(cirq._compat.proper_repr(w) for w in self.weights)) + 609 ('' if self.exponent == 1 else 610 (', exponent=' + cirq._compat.proper_repr(self.exponent))) + 611 ('' if self._global_shift == 0 else 612 (', global_shift=' + 613 cirq._compat.proper_repr(self._global_shift))) + ')') 614 615 @property 616 def qubit_generator_matrix(self): 617 generator = np.zeros((8, 8), dtype=np.complex128) 618 # w0 |110><101| + h.c. 619 generator[6, 5] = self.weights[0] 620 generator[5, 6] = self.weights[0].conjugate() 621 # w1 |110><011| + h.c. 622 generator[6, 3] = self.weights[1] 623 generator[3, 6] = self.weights[1].conjugate() 624 # w2 |101><011| + h.c. 625 generator[5, 3] = self.weights[2] 626 generator[3, 5] = self.weights[2].conjugate() 627 return generator 628 629 @staticmethod 630 def fermion_generator_components(): 631 return ( 632 ops.FermionOperator(((0, 1), (0, 0), (1, 1), (2, 0))), 633 ops.FermionOperator(((0, 1), (1, 1), (1, 0), (2, 0)), -1), 634 ops.FermionOperator(((0, 1), (1, 0), (2, 1), (2, 0))), 635 ) 636 637 @classmethod 638 def from_interaction_operator( 639 cls, 640 *, 641 operator: 'openfermion.InteractionOperator', 642 modes: Optional[Sequence[int]] = None, 643 ) -> Optional['CubicFermionicSimulationGate']: 644 if modes is None: 645 modes = (0, 1, 2) 646 i, j, k = modes 647 weights = tuple( 648 sgn * (operator.two_body_tensor[p, q, p, r] - 649 operator.two_body_tensor[p, q, r, p] - 650 operator.two_body_tensor[q, p, p, r] + 651 operator.two_body_tensor[q, p, r, p]) 652 for sgn, (p, q, 653 r) in zip([1, -1, 1], [(i, j, k), (j, i, k), (k, i, j)])) 654 if any(weights): 655 return cls(cast(Tuple[complex, complex, complex], weights)) 656 return None 657 658 def interaction_operator_generator( 659 self, 660 *, 661 operator: Optional['openfermion.InteractionOperator'] = None, 662 modes: Optional[Sequence[int]] = None 663 ) -> 'openfermion.InteractionOperator': 664 if modes is None: 665 modes = (0, 1, 2) 666 if operator is None: 667 n_modes = max(modes) + 1 668 operator = ops.InteractionOperator.zero(n_modes) 669 670 weights = tuple(w * self._exponent for w in self.weights) 671 operator.constant += self._exponent * self._global_shift 672 p, q, r = modes 673 operator.two_body_tensor[p, q, p, r] += weights[0] 674 operator.two_body_tensor[p, r, p, q] += weights[0].conjugate() 675 operator.two_body_tensor[p, q, q, r] += weights[1] 676 operator.two_body_tensor[q, r, p, q] += weights[1].conjugate() 677 operator.two_body_tensor[p, r, q, r] += weights[2] 678 operator.two_body_tensor[q, r, p, r] += weights[2].conjugate() 679 return operator 680 681 def fswap(self, i: int): 682 if i == 0: 683 self.weights = (-self.weights[1], -self.weights[0], 684 self.weights[2].conjugate()) 685 elif i == 1: 686 self.weights = (self.weights[0].conjugate(), -self.weights[2], 687 -self.weights[1]) 688 else: 689 raise ValueError(f'{i} not in (0, 1)') 690 691 692class QuarticFermionicSimulationGate(InteractionOperatorFermionicGate, 693 cirq.EigenGate): 694 r"""Rotates Hamming-weight 2 states into their bitwise complements. 695 696 With weights $(w_0, w_1, w_2)$ and exponent $t$, this gate's 697 matrix is defined as 698 699 $$ 700 e^{-i t H}, 701 $$ 702 703 where 704 705 $$ 706 H = \left(w_0 \left| 1001 \right\rangle\left\langle 0110 \right| + 707 \text{h.c.}\right) + 708 \left(w_1 \left| 1010 \right\rangle\left\langle 0101 \right| + 709 \text{h.c.}\right) + 710 \left(w_2 \left| 1100 \right\rangle\left\langle 0011 \right| + 711 \text{h.c.}\right) 712 $$ 713 714 This corresponds to the Jordan-Wigner transform of 715 716 $$ 717 H = -\left(w_0 a^{\dagger}_i a^{\dagger}_{i+3} a_{i+1} a_{i+2} + 718 \text{h.c.}\right) - 719 \left(w_1 a^{\dagger}_i a^{\dagger}_{i+2} a_{i+1} a_{i+3} + 720 \text{h.c.}\right) - 721 \left(w_2 a^{\dagger}_i a^{\dagger}_{i+1} a_{i+2} a_{i+3} + 722 \text{h.c.}\right), 723 $$ 724 725 where $a_i$, ..., $a_{i+3}$ are the annihilation operators for 726 the fermionic modes $i$, ..., $(i+3)$, respectively 727 mapped to the four qubits on which this gate acts. 728 729 730 Args: 731 weights: The weights of the terms in the Hamiltonian. 732 """ 733 734 @classmethod 735 def num_weights(cls): 736 return 3 737 738 def num_qubits(self): 739 return 4 740 741 def _eigen_components(self): 742 # projector onto subspace spanned by basis states with 743 # Hamming weight != 2 744 zero_component = np.diag( 745 [int(bin(i).count('1') != 2) for i in range(16)]) 746 747 state_pairs = (('0110', '1001'), ('0101', '1010'), ('0011', '1100')) 748 749 plus_minus_components = tuple( 750 (-abs(weight) * sign / np.pi, 751 state_swap_eigen_component(state_pair[0], state_pair[1], sign, 752 np.angle(weight))) 753 for weight, state_pair in zip(self.weights, state_pairs) 754 for sign in (-1, 1)) 755 756 return ((0, zero_component),) + plus_minus_components 757 758 def _with_exponent(self, exponent: Union[sympy.Symbol, float] 759 ) -> 'QuarticFermionicSimulationGate': 760 gate = QuarticFermionicSimulationGate(self.weights) 761 gate._exponent = exponent 762 return gate 763 764 def _decompose_(self, qubits): 765 """The goal is to effect a rotation around an axis in the XY plane in 766 each of three orthogonal 2-dimensional subspaces. 767 768 First, the following basis change is performed: 769 0000 ↦ 0001 0001 ↦ 1111 770 1111 ↦ 0010 1110 ↦ 1100 771 0010 ↦ 0000 772 0110 ↦ 0101 1101 ↦ 0011 773 1001 ↦ 0110 0100 ↦ 0100 774 1010 ↦ 1001 1011 ↦ 0111 775 0101 ↦ 1010 1000 ↦ 1000 776 1100 ↦ 1101 0111 ↦ 1011 777 0011 ↦ 1110 778 779 Note that for each 2-dimensional subspace of interest, the first two 780 qubits are the same and the right two qubits are different. The desired 781 rotations thus can be effected by a complex-version of a partial SWAP 782 gate on the latter two qubits, controlled on the first two qubits. This 783 partial SWAP-like gate can be decomposed such that it is parameterized 784 solely by a rotation in the ZY plane on the third qubit. These are the 785 `individual_rotations`; call them U0, U1, U2. 786 787 To decompose the double controlled rotations, we use four other 788 rotations V0, V1, V2, V3 (the `combined_rotations`) such that 789 U0 = V3 · V1 · V0 790 U1 = V3 · V2 · V1 791 U2 = V2 · V0 792 """ 793 794 if self._is_parameterized_(): 795 return NotImplemented 796 797 individual_rotations = [ 798 la.expm(0.5j * self.exponent * 799 np.array([[np.real(w), 1j * s * np.imag(w)], 800 [-1j * s * np.imag(w), -np.real(w)]])) 801 for s, w in zip([1, -1, -1], self.weights) 802 ] 803 804 combined_rotations = {} 805 combined_rotations[0] = la.sqrtm( 806 np.linalg.multi_dot([ 807 la.inv(individual_rotations[1]), individual_rotations[0], 808 individual_rotations[2] 809 ])) 810 combined_rotations[1] = la.inv(combined_rotations[0]) 811 combined_rotations[2] = np.linalg.multi_dot([ 812 la.inv(individual_rotations[0]), individual_rotations[1], 813 combined_rotations[0] 814 ]) 815 combined_rotations[3] = individual_rotations[0] 816 817 controlled_rotations = { 818 i: cirq.ControlledGate( 819 cirq.MatrixGate(combined_rotations[i], qid_shape=(2,))) 820 for i in range(4) 821 } 822 823 a, b, c, d = qubits 824 825 basis_change = list( 826 cirq.flatten_op_tree([ 827 cirq.CNOT(b, a), 828 cirq.CNOT(c, b), 829 cirq.CNOT(d, c), 830 cirq.CNOT(c, b), 831 cirq.CNOT(b, a), 832 cirq.CNOT(a, b), 833 cirq.CNOT(b, c), 834 cirq.CNOT(a, b), 835 [cirq.X(c), cirq.X(d)], 836 [cirq.CNOT(c, d), cirq.CNOT(d, c)], 837 [cirq.X(c), cirq.X(d)], 838 ])) 839 840 controlled_rotations = list( 841 cirq.flatten_op_tree([ 842 controlled_rotations[0](b, c), 843 cirq.CNOT(a, b), controlled_rotations[1](b, c), 844 cirq.CNOT(b, a), 845 cirq.CNOT(a, b), controlled_rotations[2](b, c), 846 cirq.CNOT(a, b), controlled_rotations[3](b, c) 847 ])) 848 849 controlled_swaps = [ 850 [cirq.CNOT(c, d), cirq.H(c)], 851 cirq.CNOT(d, c), 852 controlled_rotations, 853 cirq.CNOT(d, c), 854 [cirq.inverse(op) for op in reversed(controlled_rotations)], 855 [cirq.H(c), cirq.CNOT(c, d)], 856 ] 857 858 return [basis_change, controlled_swaps, basis_change[::-1]] 859 860 @classmethod 861 def wire_symbol(cls, use_unicode: bool): 862 return '⇊⇈' if use_unicode else 'a*a*aa' 863 864 def _apply_unitary_(self, 865 args: cirq.ApplyUnitaryArgs) -> Optional[np.ndarray]: 866 if cirq.is_parameterized(self): 867 return NotImplemented 868 869 am, bm, cm = (la.expm(-1j * self.exponent * 870 np.array([[0, w], [w.conjugate(), 0]])) 871 for w in self.weights) 872 873 a1 = args.subspace_index(0b1001) 874 b1 = args.subspace_index(0b0101) 875 c1 = args.subspace_index(0b0011) 876 877 a2 = args.subspace_index(0b0110) 878 b2 = args.subspace_index(0b1010) 879 c2 = args.subspace_index(0b1100) 880 881 cirq.apply_matrix_to_slices(args.target_tensor, 882 am, 883 slices=[a1, a2], 884 out=args.available_buffer) 885 cirq.apply_matrix_to_slices(args.available_buffer, 886 bm, 887 slices=[b1, b2], 888 out=args.target_tensor) 889 return cirq.apply_matrix_to_slices(args.target_tensor, 890 cm, 891 slices=[c1, c2], 892 out=args.available_buffer) 893 894 def __repr__(self): 895 return ('openfermion.QuarticFermionicSimulationGate(({}), ' 896 'absorb_exponent=False, ' 897 'exponent={})'.format( 898 ', '.join( 899 cirq._compat.proper_repr(v) for v in self.weights), 900 cirq._compat.proper_repr(self.exponent))) 901 902 @property 903 def qubit_generator_matrix(self): 904 """The (Hermitian) matrix G such that the gate's unitary is 905 exp(-i * G). 906 """ 907 908 generator = np.zeros((1 << 4,) * 2, dtype=np.complex128) 909 910 # w0 |1001><0110| + h.c. 911 generator[9, 6] = self.weights[0] 912 generator[6, 9] = self.weights[0].conjugate() 913 # w1 |1010><0101| + h.c. 914 generator[10, 5] = self.weights[1] 915 generator[5, 10] = self.weights[1].conjugate() 916 # w2 |1100><0011| + h.c. 917 generator[12, 3] = self.weights[2] 918 generator[3, 12] = self.weights[2].conjugate() 919 return generator 920 921 @staticmethod 922 def fermion_generator_components(): 923 return ( 924 ops.FermionOperator(((0, 1), (1, 0), (2, 0), (3, 1)), -1), 925 ops.FermionOperator(((0, 1), (1, 0), (2, 1), (3, 0))), 926 ops.FermionOperator(((0, 1), (1, 1), (2, 0), (3, 0)), -1), 927 ) 928 929 @classmethod 930 def from_interaction_operator( 931 cls, 932 *, 933 operator: 'openfermion.InteractionOperator', 934 modes: Optional[Sequence[int]] = None, 935 ) -> Optional['QuarticFermionicSimulationGate']: 936 if modes is None: 937 modes = (0, 1, 2, 3) 938 i, j, k, l = modes 939 weights = tuple( 940 (operator.two_body_tensor[p, q, r, s] - 941 operator.two_body_tensor[p, q, s, r] - 942 operator.two_body_tensor[q, p, r, s] + 943 operator.two_body_tensor[q, p, s, r]) 944 for p, q, r, s in [(i, l, j, k), (i, k, j, l), (i, j, k, l)]) 945 if any(weights): 946 return cls(cast(Tuple[complex, complex, complex], weights)) 947 return None 948 949 def interaction_operator_generator( 950 self, 951 operator: Optional['openfermion.InteractionOperator'] = None, 952 modes: Optional[Sequence[int]] = None 953 ) -> 'openfermion.InteractionOperator': 954 if modes is None: 955 modes = (0, 1, 2, 3) 956 if operator is None: 957 n_modes = max(modes) + 1 958 operator = ops.InteractionOperator.zero(n_modes) 959 960 weights = tuple(w * self._exponent for w in self.weights) 961 operator.constant += self._exponent * self._global_shift 962 p, q, r, s = modes 963 operator.two_body_tensor[p, s, q, r] += weights[0] 964 operator.two_body_tensor[q, r, p, s] += weights[0].conjugate() 965 operator.two_body_tensor[p, r, q, s] += weights[1] 966 operator.two_body_tensor[q, s, p, r] += weights[1].conjugate() 967 operator.two_body_tensor[p, q, r, s] += weights[2] 968 operator.two_body_tensor[r, s, p, q] += weights[2].conjugate() 969 return operator 970 971 def fswap(self, i: int): 972 if i == 0: 973 self.weights = (self.weights[1].conjugate(), 974 self.weights[0].conjugate(), -self.weights[2]) 975 elif i == 1: 976 self.weights = (-self.weights[0], self.weights[2], self.weights[1]) 977 elif i == 2: 978 self.weights = (self.weights[1], self.weights[0], -self.weights[2]) 979 else: 980 raise ValueError(f'{i} not in (0, 1, 2)') 981