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"""Bravyi-Kitaev transform on fermionic operators.""" 13 14from openfermion.ops.operators import (FermionOperator, MajoranaOperator, 15 QubitOperator) 16from openfermion.ops.representations import InteractionOperator 17from openfermion.utils.operator_utils import count_qubits 18 19 20def bravyi_kitaev(operator, n_qubits=None): 21 """Apply the Bravyi-Kitaev transform. 22 23 Implementation from arXiv:quant-ph/0003137 and 24 "A New Data Structure for Cumulative Frequency Tables" by Peter M. Fenwick. 25 26 Note that this implementation is equivalent to the one described in 27 arXiv:1208.5986, and is different from the one described in 28 arXiv:1701.07072. The one described in arXiv:1701.07072 is implemented 29 in OpenFermion as `bravyi_kitaev_tree`. 30 31 Args: 32 operator (openfermion.ops.FermionOperator): 33 A FermionOperator to transform. 34 n_qubits (int|None): 35 Can force the number of qubits in the resulting operator above the 36 number that appear in the input operator. 37 38 Returns: 39 transformed_operator: An instance of the QubitOperator class. 40 41 Raises: 42 ValueError: Invalid number of qubits specified. 43 """ 44 if isinstance(operator, FermionOperator): 45 return _bravyi_kitaev_fermion_operator(operator, n_qubits) 46 if isinstance(operator, MajoranaOperator): 47 return _bravyi_kitaev_majorana_operator(operator, n_qubits) 48 if isinstance(operator, InteractionOperator): 49 return _bravyi_kitaev_interaction_operator(operator, n_qubits) 50 raise TypeError("Couldn't apply the Bravyi-Kitaev Transform to object " 51 "of type {}.".format(type(operator))) 52 53 54def _update_set(index, n_qubits): 55 """The bits that need to be updated upon flipping the occupancy 56 of a mode.""" 57 indices = set() 58 59 # For bit manipulation we need to count from 1 rather than 0 60 index += 1 61 # Ensure index is not a member of the set 62 index += index & -index 63 64 while index <= n_qubits: 65 indices.add(index - 1) 66 # Add least significant one to index 67 # E.g. 00010100 -> 00011000 68 index += index & -index 69 return indices 70 71 72def _occupation_set(index): 73 """The bits whose parity stores the occupation of mode `index`.""" 74 indices = set() 75 76 # For bit manipulation we need to count from 1 rather than 0 77 index += 1 78 79 indices.add(index - 1) 80 parent = index & (index - 1) 81 index -= 1 82 while index != parent: 83 indices.add(index - 1) 84 # Remove least significant one from index 85 # E.g. 00010100 -> 00010000 86 index &= index - 1 87 return indices 88 89 90def _parity_set(index): 91 """The bits whose parity stores the parity of the bits 0 .. `index`.""" 92 indices = set() 93 94 while index > 0: 95 indices.add(index - 1) 96 # Remove least significant one from index 97 # E.g. 00010100 -> 00010000 98 index &= index - 1 99 return indices 100 101 102def _bravyi_kitaev_majorana_operator(operator, n_qubits): 103 # Compute the number of qubits. 104 N = count_qubits(operator) 105 if n_qubits is None: 106 n_qubits = N 107 if n_qubits < N: 108 raise ValueError('Invalid number of qubits specified.') 109 110 # Compute transformed operator. 111 transformed_terms = (_transform_majorana_term(term=term, 112 coefficient=coeff, 113 n_qubits=n_qubits) 114 for term, coeff in operator.terms.items()) 115 return inline_sum(summands=transformed_terms, seed=QubitOperator()) 116 117 118def _transform_majorana_term(term, coefficient, n_qubits): 119 # Build the Bravyi-Kitaev transformed operators. 120 transformed_ops = (_transform_majorana_operator(majorana_index, n_qubits) 121 for majorana_index in term) 122 return inline_product(factors=transformed_ops, 123 seed=QubitOperator((), coefficient)) 124 125 126def _transform_majorana_operator(majorana_index, n_qubits): 127 q, b = divmod(majorana_index, 2) 128 129 update_set = _update_set(q, n_qubits) 130 update_set.add(q) 131 occupation_set = _occupation_set(q) 132 parity_set = _parity_set(q) 133 134 if b: 135 return QubitOperator([(q, 'Y')] + [(i, 'X') for i in update_set - {q}] + 136 [(i, 'Z') 137 for i in (parity_set ^ occupation_set) - {q}]) 138 else: 139 return QubitOperator([(i, 'X') for i in update_set] + 140 [(i, 'Z') for i in parity_set]) 141 142 143def _transform_operator_term(term, coefficient, n_qubits): 144 """ 145 Args: 146 term (list[tuple[int, int]]): 147 A list of (mode, raising-vs-lowering) ladder operator terms. 148 coefficient (float): 149 n_qubits (int): 150 Returns: 151 QubitOperator: 152 """ 153 154 # Build the Bravyi-Kitaev transformed operators. 155 transformed_ladder_ops = (_transform_ladder_operator( 156 ladder_operator, n_qubits) for ladder_operator in term) 157 return inline_product(factors=transformed_ladder_ops, 158 seed=QubitOperator((), coefficient)) 159 160 161def _bravyi_kitaev_fermion_operator(operator, n_qubits): 162 # Compute the number of qubits. 163 N = count_qubits(operator) 164 if n_qubits is None: 165 n_qubits = N 166 if n_qubits < N: 167 raise ValueError('Invalid number of qubits specified.') 168 169 # Compute transformed operator. 170 transformed_terms = (_transform_operator_term( 171 term=term, coefficient=operator.terms[term], n_qubits=n_qubits) 172 for term in operator.terms) 173 return inline_sum(summands=transformed_terms, seed=QubitOperator()) 174 175 176def _transform_ladder_operator(ladder_operator, n_qubits): 177 """ 178 Args: 179 ladder_operator (tuple[int, int]): the ladder operator 180 n_qubits (int): the number of qubits 181 Returns: 182 QubitOperator 183 """ 184 index, action = ladder_operator 185 186 update_set = _update_set(index, n_qubits) 187 update_set.add(index) 188 occupation_set = _occupation_set(index) 189 parity_set = _parity_set(index) 190 191 # Initialize the transformed majorana operator (a_p^\dagger + a_p) / 2 192 transformed_operator = QubitOperator([(i, 'X') for i in update_set] + 193 [(i, 'Z') for i in parity_set], .5) 194 # Get the transformed (a_p^\dagger - a_p) / 2 195 # Below is equivalent to X(update_set) * Z(parity_set ^ occupation_set) 196 transformed_majorana_difference = QubitOperator( 197 [(index, 'Y')] + [(i, 'X') for i in update_set - {index}] + 198 [(i, 'Z') for i in (parity_set ^ occupation_set) - {index}], -.5j) 199 200 # Raising 201 if action == 1: 202 transformed_operator += transformed_majorana_difference 203 # Lowering 204 else: 205 transformed_operator -= transformed_majorana_difference 206 207 return transformed_operator 208 209 210def inline_sum(summands, seed): 211 """Computes a sum, using the __iadd__ operator. 212 Args: 213 seed (T): The starting total. The zero value. 214 summands (iterable[T]): Values to add (with +=) into the total. 215 Returns: 216 T: The result of adding all the factors into the zero value. 217 """ 218 for r in summands: 219 seed += r 220 return seed 221 222 223def inline_product(factors, seed): 224 """Computes a product, using the __imul__ operator. 225 Args: 226 seed (T): The starting total. The unit value. 227 factors (iterable[T]): Values to multiply (with *=) into the total. 228 Returns: 229 T: The result of multiplying all the factors into the unit value. 230 """ 231 for r in factors: 232 seed *= r 233 return seed 234 235 236def _bravyi_kitaev_interaction_operator(interaction_operator, n_qubits): 237 """Implementation of the Bravyi-Kitaev transformation for OpenFermion 238 Interaction Operators. This implementation is equivalent to that described 239 in arXiv:1208.5986, and has been written to optimize compute time by using 240 algebraic expressions for general products a_i^dagger a_j^dagger as outlined 241 in Table II of Seeley, Richard, Love. 242 """ 243 244 one_body = interaction_operator.one_body_tensor 245 two_body = interaction_operator.two_body_tensor 246 constant_term = interaction_operator.constant 247 248 # Compute the number of qubits. 249 N = len(one_body) 250 if n_qubits is None: 251 n_qubits = N 252 if n_qubits < N: 253 raise ValueError('Invalid number of qubits specified.') 254 255 qubit_hamiltonian = QubitOperator() 256 qubit_hamiltonian_op = [] 257 qubit_hamiltonian_coef = [] 258 259 # For cases A - F see Table II of Seeley, Richard, Love 260 for i in range(n_qubits): 261 # A. Number operators: n_i 262 if abs(one_body[i, i]) > 0: 263 qubit_hamiltonian += _qubit_operator_creation( 264 *_seeley_richard_love(i, i, one_body[i, i], n_qubits)) 265 266 for j in range(i): 267 # Case B: Coulomb and exchange operators 268 if abs(one_body[i, j]) > 0: 269 operators, coef_list = _seeley_richard_love( 270 i, j, one_body[i, j], n_qubits) 271 qubit_hamiltonian_op.extend(operators) 272 qubit_hamiltonian_coef.extend(coef_list) 273 274 operators, coef_list = _seeley_richard_love( 275 j, i, one_body[i, j].conj(), n_qubits) 276 qubit_hamiltonian_op.extend(operators) 277 qubit_hamiltonian_coef.extend(coef_list) 278 279 coef = _two_body_coef(two_body, i, j, j, i) / 4 280 if abs(coef) > 0: 281 qubit_hamiltonian_op.append( 282 tuple((index, "Z") for index in _occupation_set(i))) 283 qubit_hamiltonian_op.append( 284 tuple((index, "Z") for index in _occupation_set(j))) 285 qubit_hamiltonian_op.append( 286 tuple((index, "Z") for index in _F_ij_set(i, j))) 287 288 qubit_hamiltonian_coef.append(-coef) 289 qubit_hamiltonian_coef.append(-coef) 290 qubit_hamiltonian_coef.append(coef) 291 constant_term += coef 292 293 # C. Number-excitation operators: n_i a_j^d a_k 294 for i in range(n_qubits): 295 for j in range(n_qubits): 296 for k in range(j): 297 if i not in (j, k): 298 coef = _two_body_coef(two_body, i, j, k, i) 299 300 if abs(coef) > 0: 301 number = _qubit_operator_creation( 302 *_seeley_richard_love(i, i, 1, n_qubits)) 303 304 excitation_op, excitation_coef = _seeley_richard_love( 305 j, k, coef, n_qubits) 306 operators_hc, coef_list_hc = _seeley_richard_love( 307 k, j, coef.conj(), n_qubits) 308 excitation_op.extend(operators_hc) 309 excitation_coef.extend(coef_list_hc) 310 excitation = _qubit_operator_creation( 311 excitation_op, excitation_coef) 312 313 number *= excitation 314 qubit_hamiltonian += number 315 316 # D. Double-excitation operators: c_i^d c_j^d c_k c_l 317 for i in range(n_qubits): 318 for j in range(i): 319 for k in range(j): 320 for l in range(k): 321 coef = -_two_body_coef(two_body, i, j, k, l) 322 if abs(coef) > 0: 323 qubit_hamiltonian += _hermitian_one_body_product( 324 i, j, k, l, coef, n_qubits) 325 326 coef = -_two_body_coef(two_body, i, k, j, l) 327 if abs(coef) > 0: 328 qubit_hamiltonian += _hermitian_one_body_product( 329 i, k, j, l, coef, n_qubits) 330 331 coef = -_two_body_coef(two_body, i, l, j, k) 332 if abs(coef) > 0: 333 qubit_hamiltonian += _hermitian_one_body_product( 334 i, l, j, k, coef, n_qubits) 335 336 qubit_hamiltonian_op.append(()) 337 qubit_hamiltonian_coef.append(constant_term) 338 qubit_hamiltonian += _qubit_operator_creation(qubit_hamiltonian_op, 339 qubit_hamiltonian_coef) 340 341 return qubit_hamiltonian 342 343 344def _two_body_coef(two_body, a, b, c, d): 345 return two_body[a, b, c, d] - two_body[a, b, d, c] + two_body[ 346 b, a, d, c] - two_body[b, a, c, d] 347 348 349def _hermitian_one_body_product(a, b, c, d, coef, n_qubits): 350 """ Takes the 4 indices for a two-body operator and constructs the 351 Bravyi-Kitaev form by splitting the two-body into 2 one-body operators, 352 multiplying them together and then re-adding the Hermitian conjugate to 353 give a Hermitian operator. """ 354 355 c_dag_c_ac = _qubit_operator_creation( 356 *_seeley_richard_love(a, c, coef, n_qubits)) 357 c_dag_c_bd = _qubit_operator_creation( 358 *_seeley_richard_love(b, d, 1, n_qubits)) 359 c_dag_c_ac *= c_dag_c_bd 360 hermitian_sum = c_dag_c_ac 361 362 c_dag_c_ca = _qubit_operator_creation( 363 *_seeley_richard_love(c, a, coef.conj(), n_qubits)) 364 c_dag_c_db = _qubit_operator_creation( 365 *_seeley_richard_love(d, b, 1, n_qubits)) 366 c_dag_c_ca *= c_dag_c_db 367 368 hermitian_sum += c_dag_c_ca 369 return hermitian_sum 370 371 372def _qubit_operator_creation(operators, coefficents): 373 """ Takes a list of tuples for operators/indices, and another for 374 coefficents""" 375 376 qubit_operator = QubitOperator() 377 378 for index in zip(operators, coefficents): 379 qubit_operator += QubitOperator(index[0], index[1]) 380 381 return qubit_operator 382 383 384def _seeley_richard_love(i, j, coef, n_qubits): 385 """Algebraic expressions for general products of the form a_i^d a_j term in 386 the Bravyi-Kitaev basis. These expressions vary in form depending on the 387 parity of the indices i and j, as well as onthe overlaps between the parity 388 and update sets of the indices""" 389 390 seeley_richard_love_op = [] 391 seeley_richard_love_coef = [] 392 coef *= 0.25 393 # Case 0 394 if i == j: # Simplifies to the number operator 395 seeley_richard_love_op.append( 396 tuple((index, "Z") for index in _occupation_set(i))) 397 seeley_richard_love_coef.append(-coef * 2) 398 399 seeley_richard_love_op.append(()) 400 seeley_richard_love_coef.append(coef * 2) 401 402 # Case 1 403 elif i % 2 == 0 and j % 2 == 0: 404 x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits)) 405 y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits)) 406 z_pad = tuple( 407 (index, "Z") for index in _P0_ij_diff_a_set(i, j, n_qubits)) 408 409 left_pad = x_pad + y_pad + z_pad 410 411 seeley_richard_love_op.append(left_pad + ((j, "Y"), (i, "X"))) 412 seeley_richard_love_op.append(left_pad + ((j, "X"), (i, "Y"))) 413 seeley_richard_love_op.append(left_pad + ((j, "X"), (i, "X"))) 414 seeley_richard_love_op.append(left_pad + ((j, "Y"), (i, "Y"))) 415 416 if i < j: 417 seeley_richard_love_coef.append(coef) 418 seeley_richard_love_coef.append(-coef) 419 seeley_richard_love_coef.append(complex(0, -coef)) 420 seeley_richard_love_coef.append(complex(0, -coef)) 421 422 else: # whenever i < j introduce phase of -j 423 seeley_richard_love_coef.append(complex(0, -coef)) 424 seeley_richard_love_coef.append(complex(0, coef)) 425 seeley_richard_love_coef.append(-coef) 426 seeley_richard_love_coef.append(-coef) 427 428 # Case 2 429 elif i % 2 == 1 and j % 2 == 0 and i not in _parity_set(j): 430 x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits)) 431 y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits)) 432 433 left_pad = x_pad + y_pad 434 435 right_pad_1 = tuple( 436 (index, "Z") 437 for index in _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits)) 438 right_pad_2 = tuple( 439 (index, "Z") 440 for index in _P2_ij_set(i, j) - _alpha_set(i, j, n_qubits)) 441 442 seeley_richard_love_op.append(left_pad + ((j, "Y"), 443 (i, "X")) + right_pad_1) 444 seeley_richard_love_op.append(left_pad + ((j, "X"), 445 (i, "X")) + right_pad_1) 446 seeley_richard_love_op.append(left_pad + ((j, "X"), 447 (i, "Y")) + right_pad_2) 448 seeley_richard_love_op.append(left_pad + ((j, "Y"), 449 (i, "Y")) + right_pad_2) 450 451 if i < j: 452 seeley_richard_love_coef.append(coef) 453 seeley_richard_love_coef.append(complex(0, -coef)) 454 seeley_richard_love_coef.append(-coef) 455 seeley_richard_love_coef.append(complex(0, -coef)) 456 457 else: 458 seeley_richard_love_coef.append(complex(0, -coef)) 459 seeley_richard_love_coef.append(-coef) 460 seeley_richard_love_coef.append(complex(0, coef)) 461 seeley_richard_love_coef.append(-coef) 462 463 # Case 3 464 elif i % 2 == 1 and j % 2 == 0 and i in _parity_set(j): 465 left_pad = tuple((index, "X") for index in _U_ij_set(i, j, n_qubits)) 466 right_pad_1 = tuple((index, "Z") for index in _P0_ij_set(i, j) - {i}) 467 right_pad_2 = tuple((index, "Z") for index in _P2_ij_set(i, j) - {i}) 468 469 seeley_richard_love_op.append(left_pad + ((j, "Y"), 470 (i, "Y")) + right_pad_1) 471 seeley_richard_love_coef.append(coef) 472 473 seeley_richard_love_op.append(left_pad + ((j, "X"), 474 (i, "Y")) + right_pad_1) 475 seeley_richard_love_coef.append(complex(0, -coef)) 476 477 seeley_richard_love_op.append(left_pad + ((j, "X"), 478 (i, "X")) + right_pad_2) 479 seeley_richard_love_coef.append(coef) 480 481 seeley_richard_love_op.append(left_pad + ((j, "Y"), 482 (i, "X")) + right_pad_2) 483 seeley_richard_love_coef.append(complex(0, coef)) 484 485 # Case 4 486 elif i % 2 == 0 and j % 2 == 1 and i not in _parity_set( 487 j) and j not in _update_set(i, n_qubits): 488 x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits)) 489 y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits)) 490 left_pad = x_pad + y_pad 491 492 right_pad_1 = tuple( 493 (index, "Z") 494 for index in _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits)) 495 right_pad_2 = tuple( 496 (index, "Z") 497 for index in _P1_ij_set(i, j) - _alpha_set(i, j, n_qubits)) 498 499 seeley_richard_love_op.append(left_pad + ((j, "X"), 500 (i, "Y")) + right_pad_1) 501 seeley_richard_love_op.append(left_pad + ((j, "X"), 502 (i, "X")) + right_pad_1) 503 seeley_richard_love_op.append(left_pad + ((j, "Y"), 504 (i, "X")) + right_pad_2) 505 seeley_richard_love_op.append(left_pad + ((j, "Y"), 506 (i, "Y")) + right_pad_2) 507 508 if i < j: 509 seeley_richard_love_coef.append(-coef) 510 seeley_richard_love_coef.append(complex(0, -coef)) 511 seeley_richard_love_coef.append(coef) 512 seeley_richard_love_coef.append(complex(0, -coef)) 513 else: 514 seeley_richard_love_coef.append(complex(0, coef)) 515 seeley_richard_love_coef.append(-coef) 516 seeley_richard_love_coef.append(complex(0, -coef)) 517 seeley_richard_love_coef.append(-coef) 518 519 # Case 5 520 elif i % 2 == 0 and j % 2 == 1 and i not in _parity_set( 521 j) and j in _update_set(i, n_qubits): 522 x_range_1 = _U_ij_set(i, j, n_qubits) - {j} 523 left_pad_1 = tuple((index, "X") for index in x_range_1) 524 525 x_range_2 = x_range_1 - _alpha_set(i, j, n_qubits) 526 left_pad_2 = tuple((index, "X") for index in x_range_2) 527 528 y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits)) 529 z_range_1 = _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits) 530 z_pad = tuple((index, "Z") for index in z_range_1) 531 right_pad_1 = y_pad + z_pad 532 533 z_range_2 = _P1_ij_set(i, j).union({j}) 534 right_pad_2 = tuple((index, "Z") for index in z_range_2) 535 536 seeley_richard_love_op.append(left_pad_2 + ((i, "Y"),) + right_pad_1) 537 seeley_richard_love_coef.append(-coef) 538 539 seeley_richard_love_op.append(left_pad_2 + ((i, "X"),) + right_pad_1) 540 seeley_richard_love_coef.append(complex( 541 0, -coef)) # Phase flip of -1 relative to original paper 542 543 seeley_richard_love_op.append(left_pad_1 + ((i, "Y"),) + right_pad_2) 544 seeley_richard_love_coef.append(complex(0, coef)) 545 546 seeley_richard_love_op.append(left_pad_1 + ((i, "X"),) + right_pad_2) 547 seeley_richard_love_coef.append(-coef) 548 549 # Case 6 550 elif i % 2 == 0 and j % 2 == 1 and i in _parity_set(j) and j in _update_set( 551 i, n_qubits): 552 left_pad = tuple( 553 (index, "X") for index in _U_ij_set(i, j, n_qubits) - {j}) 554 right_pad = tuple((index, "Z") for index in _P1_ij_set(i, j).union({j})) 555 556 seeley_richard_love_op.append(left_pad + ((i, "X"),)) 557 seeley_richard_love_coef.append(coef) 558 559 seeley_richard_love_op.append(left_pad + ((i, "Y"),)) 560 seeley_richard_love_coef.append(complex(0, -coef)) 561 562 seeley_richard_love_op.append(left_pad + ((i, "Y"),) + right_pad) 563 seeley_richard_love_coef.append(complex(0, coef)) 564 565 seeley_richard_love_op.append(left_pad + ((i, "X"),) + right_pad) 566 seeley_richard_love_coef.append(-coef) 567 568 # Case 7 569 elif i % 2 == 1 and j % 2 == 1 and i not in _parity_set( 570 j) and j not in _update_set(i, n_qubits): 571 x_pad = tuple((index, "X") for index in _U_diff_a_set(i, j, n_qubits)) 572 y_pad = tuple((index, "Y") for index in _alpha_set(i, j, n_qubits)) 573 left_pad = x_pad + y_pad 574 575 z_range_1 = _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits) 576 right_pad_1 = tuple((index, "Z") for index in z_range_1) 577 578 z_range_2 = _P1_ij_set(i, j) - _alpha_set(i, j, n_qubits) 579 right_pad_2 = tuple((index, "Z") for index in z_range_2) 580 581 z_range_3 = _P2_ij_set(i, j) - _alpha_set(i, j, n_qubits) 582 right_pad_3 = tuple((index, "Z") for index in z_range_3) 583 584 z_range_4 = _P3_ij_set(i, j) - _alpha_set(i, j, n_qubits) 585 right_pad_4 = tuple((index, "Z") for index in z_range_4) 586 587 if i < j: 588 seeley_richard_love_op.append(left_pad + ((j, "X"), 589 (i, "X")) + right_pad_1) 590 seeley_richard_love_coef.append(complex(0, -coef)) 591 592 seeley_richard_love_op.append(left_pad + ((j, "Y"), 593 (i, "X")) + right_pad_2) 594 seeley_richard_love_coef.append(coef) 595 596 seeley_richard_love_op.append(left_pad + ((j, "X"), 597 (i, "Y")) + right_pad_3) 598 seeley_richard_love_coef.append(-coef) 599 600 seeley_richard_love_op.append(left_pad + ((j, "Y"), 601 (i, "Y")) + right_pad_4) 602 seeley_richard_love_coef.append(complex(0, -coef)) 603 604 else: 605 seeley_richard_love_op.append(left_pad + ((j, "X"), 606 (i, "X")) + right_pad_1) 607 seeley_richard_love_coef.append(-coef) 608 609 seeley_richard_love_op.append(left_pad + ((j, "Y"), 610 (i, "X")) + right_pad_2) 611 seeley_richard_love_coef.append(complex(0, -coef)) 612 613 seeley_richard_love_op.append(left_pad + ((j, "X"), 614 (i, "Y")) + right_pad_3) 615 seeley_richard_love_coef.append(complex(0, coef)) 616 617 seeley_richard_love_op.append(left_pad + ((j, "Y"), 618 (i, "Y")) + right_pad_4) 619 seeley_richard_love_coef.append(-coef) 620 621 # Case 8 622 elif i % 2 == 1 and j % 2 == 1 and i in _parity_set( 623 j) and j not in _update_set(i, n_qubits): 624 left_pad = tuple((index, "X") for index in _U_ij_set(i, j, n_qubits)) 625 626 z_range_1 = _P0_ij_set(i, j) - {i} 627 right_pad_1 = tuple((index, "Z") for index in z_range_1) 628 629 z_range_2 = _P1_ij_set(i, j) - {i} 630 right_pad_2 = tuple((index, "Z") for index in z_range_2) 631 632 z_range_3 = _P2_ij_set(i, j) - {i} 633 right_pad_3 = tuple((index, "Z") for index in z_range_3) 634 635 z_range_4 = _P3_ij_set(i, j) - {i} 636 right_pad_4 = tuple((index, "Z") for index in z_range_4) 637 638 seeley_richard_love_op.append(left_pad + ((j, "X"), 639 (i, "Y")) + right_pad_1) 640 seeley_richard_love_coef.append(complex(0, -coef)) 641 642 seeley_richard_love_op.append(left_pad + ((j, "Y"), 643 (i, "Y")) + right_pad_2) 644 seeley_richard_love_coef.append(coef) 645 646 seeley_richard_love_op.append(left_pad + ((j, "X"), 647 (i, "X")) + right_pad_3) 648 seeley_richard_love_coef.append(coef) 649 650 seeley_richard_love_op.append(left_pad + ((j, "Y"), 651 (i, "X")) + right_pad_4) 652 seeley_richard_love_coef.append(complex(0, coef)) 653 654 # Case 9 655 elif i % 2 == 1 and j % 2 == 1 and i not in _parity_set( 656 j) and j in _update_set(i, n_qubits): 657 x_range_1 = _U_ij_set(i, j, n_qubits) - {j} 658 left_pad_3 = tuple((index, "X") for index in x_range_1) 659 660 x_range_2 = x_range_1 - _alpha_set(i, j, n_qubits) 661 left_pad_1 = tuple((index, "X") for index in x_range_2) 662 663 x_range_3 = x_range_1.union({i}) - _alpha_set(i, j, n_qubits) 664 left_pad_2 = tuple((index, "X") for index in x_range_3) 665 666 z_range_1 = _P2_ij_set(i, j) - _alpha_set(i, j, n_qubits) 667 z_pad_1 = tuple((index, "Z") for index in z_range_1) 668 y_range = _alpha_set(i, j, n_qubits) 669 y_pad = tuple((index, "Y") for index in y_range) 670 right_pad_1 = z_pad_1 + y_pad 671 672 z_range_2 = _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits) 673 z_pad_2 = tuple((index, "Z") for index in z_range_2) 674 right_pad_2 = z_pad_2 + y_pad 675 676 z_range_3 = _P1_ij_set(i, j).union({j}) 677 right_pad_3 = tuple((index, "Z") for index in z_range_3) 678 679 z_range_4 = _P3_ij_set(i, j).union({j}) 680 right_pad_4 = tuple((index, "Z") for index in z_range_4) 681 682 seeley_richard_love_op.append(left_pad_1 + ((i, "Y"),) + right_pad_1) 683 seeley_richard_love_coef.append( 684 -coef) # phase of -j relative to original paper 685 686 seeley_richard_love_op.append(left_pad_2 + right_pad_2) 687 seeley_richard_love_coef.append(complex( 688 0, -coef)) # phase of -j relative to original paper 689 690 seeley_richard_love_op.append(left_pad_3 + ((i, "X"),) + right_pad_3) 691 seeley_richard_love_coef.append(-coef) 692 693 seeley_richard_love_op.append(left_pad_3 + ((i, "Y"),) + right_pad_4) 694 seeley_richard_love_coef.append(complex(0, coef)) 695 696 # Case 10 697 elif i % 2 == 1 and j % 2 == 1 and i in _parity_set(j) and j in _update_set( 698 i, n_qubits): 699 left_pad = tuple( 700 (index, "X") for index in _U_ij_set(i, j, n_qubits) - {j}) 701 right_pad_1 = tuple((index, "Z") for index in _P0_ij_set(i, j) - {i}) 702 right_pad_2 = tuple((index, "Z") for index in _P2_ij_set(i, j) - {i}) 703 right_pad_3 = tuple((index, "Z") for index in _P1_ij_set(i, j)) 704 right_pad_4 = tuple((index, "Z") for index in _P3_ij_set(i, j)) 705 706 seeley_richard_love_op.append(left_pad + ((i, "Y"),) + right_pad_1) 707 seeley_richard_love_coef.append(complex(0, -coef)) 708 709 seeley_richard_love_op.append(left_pad + ((i, "X"),) + right_pad_2) 710 seeley_richard_love_coef.append(coef) 711 712 seeley_richard_love_op.append(left_pad + ((j, "Z"), 713 (i, "X")) + right_pad_3) 714 seeley_richard_love_coef.append(-coef) 715 716 seeley_richard_love_op.append(left_pad + ((j, "Z"), 717 (i, "Y")) + right_pad_4) 718 seeley_richard_love_coef.append(complex(0, coef)) 719 720 return seeley_richard_love_op, seeley_richard_love_coef 721 722 723def _remainder_set(index): 724 return _parity_set(index) - _occupation_set(index) 725 726 727def _F_ij_set(i, j): 728 return _occupation_set(i).symmetric_difference(_occupation_set(j)) 729 730 731def _P0_ij_set(i, j): 732 """ The symmetric difference of sets P(i) and P(j). """ 733 return _parity_set(i).symmetric_difference(_parity_set(j)) 734 735 736def _P1_ij_set(i, j): 737 """ The symmetric difference of sets P(i) and R(j). """ 738 return _parity_set(i).symmetric_difference(_remainder_set(j)) 739 740 741def _P2_ij_set(i, j): 742 """ The symmetric difference of sets R(i) and P(j). """ 743 return _remainder_set(i).symmetric_difference(_parity_set(j)) 744 745 746def _P3_ij_set(i, j): 747 """ The symmetric difference of sets R(i) and R(j). """ 748 return _remainder_set(i).symmetric_difference(_remainder_set(j)) 749 750 751def _U_ij_set(i, j, n_qubits): 752 """ The symmetric difference of sets U(i) and U(j)""" 753 return _update_set(i, 754 n_qubits).symmetric_difference(_update_set(j, n_qubits)) 755 756 757def _U_diff_a_set(i, j, n_qubits): 758 """ Calculates the set {member U_ij diff alpha_ij}. """ 759 return _U_ij_set(i, j, n_qubits) - _alpha_set(i, j, n_qubits) 760 761 762def _P0_ij_diff_a_set(i, j, n_qubits): 763 """ Calculates the set {member P_ij diff alpha_ij}. """ 764 return _P0_ij_set(i, j) - _alpha_set(i, j, n_qubits) 765 766 767def _alpha_set(i, j, n_qubits): 768 return _update_set(i, n_qubits).intersection(_parity_set(j)) 769