From 5503237c751685773b0f514d5a059a30c29eb102 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Thu, 19 Dec 2024 12:57:53 -0600 Subject: [PATCH] Bravyi-Kitaev implementation (#35) This is the implementation of the Bravyi-Kitaev fermionic transformation. The was largely a translation from the OpenFermion BK python implementation into C++: https://github.com/quantumlib/OpenFermion/blob/master/src/openfermion/transforms/opconversions/bravyi_kitaev.py --- .../fermion_compilers/bravyi_kitaev.h | 29 + libs/solvers/lib/CMakeLists.txt | 1 + .../fermion_compilers/bravyi_kitaev.cpp | 699 ++++++++++++++++++ .../python/bindings/solvers/py_solvers.cpp | 152 ++++ libs/solvers/python/tests/test_molecule.py | 47 ++ libs/solvers/unittests/CMakeLists.txt | 7 +- libs/solvers/unittests/support/h2_pyscf_hf.py | 73 ++ libs/solvers/unittests/test_bravyi_kitaev.cpp | 243 ++++++ 8 files changed, 1250 insertions(+), 1 deletion(-) create mode 100644 libs/solvers/include/cudaq/solvers/operators/molecule/fermion_compilers/bravyi_kitaev.h create mode 100644 libs/solvers/lib/operators/molecule/fermion_compilers/bravyi_kitaev.cpp create mode 100644 libs/solvers/unittests/support/h2_pyscf_hf.py create mode 100644 libs/solvers/unittests/test_bravyi_kitaev.cpp diff --git a/libs/solvers/include/cudaq/solvers/operators/molecule/fermion_compilers/bravyi_kitaev.h b/libs/solvers/include/cudaq/solvers/operators/molecule/fermion_compilers/bravyi_kitaev.h new file mode 100644 index 0000000..a999647 --- /dev/null +++ b/libs/solvers/include/cudaq/solvers/operators/molecule/fermion_compilers/bravyi_kitaev.h @@ -0,0 +1,29 @@ +/******************************************************************************* + * Copyright (c) 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ +#pragma once + +#include "cudaq/solvers/operators/molecule/fermion_compiler.h" + +namespace cudaq::solvers { + +/// @brief Helper function used by the Bravyi-Kitaev transformation. +cudaq::spin_op seeley_richard_love(std::size_t i, std::size_t j, + std::complex coef, int n_qubits); + +/// @brief Map fermionic operators to spin operators via the +/// Bravyi-Kitaev transformation. +class bravyi_kitaev : public fermion_compiler { +public: + cudaq::spin_op generate(const double constant, const cudaqx::tensor<> &hpq, + const cudaqx::tensor<> &hpqrs, + const cudaqx::heterogeneous_map &options) override; + + CUDAQ_EXTENSION_CREATOR_FUNCTION(fermion_compiler, bravyi_kitaev) +}; +CUDAQ_REGISTER_TYPE(bravyi_kitaev) +} // namespace cudaq::solvers diff --git a/libs/solvers/lib/CMakeLists.txt b/libs/solvers/lib/CMakeLists.txt index 2fcb44e..019b4ac 100644 --- a/libs/solvers/lib/CMakeLists.txt +++ b/libs/solvers/lib/CMakeLists.txt @@ -15,6 +15,7 @@ add_library(cudaq-solvers SHARED operators/molecule/drivers/pyscf_driver.cpp operators/molecule/fermion_compilers/fermion_compiler.cpp operators/molecule/fermion_compilers/jordan_wigner.cpp + operators/molecule/fermion_compilers/bravyi_kitaev.cpp operators/molecule/molecule.cpp operators/graph/max_cut.cpp operators/graph/clique.cpp diff --git a/libs/solvers/lib/operators/molecule/fermion_compilers/bravyi_kitaev.cpp b/libs/solvers/lib/operators/molecule/fermion_compilers/bravyi_kitaev.cpp new file mode 100644 index 0000000..02f3cce --- /dev/null +++ b/libs/solvers/lib/operators/molecule/fermion_compilers/bravyi_kitaev.cpp @@ -0,0 +1,699 @@ +/******************************************************************************* + * Copyright (c) 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +/****************************************************************************** + * * + * This file was translated and modified from bravyi_kitaev.py * + * which was adapted from https://doi.org/10.1063/1.4768229 * + * Original work Copyright OpenFermion * + * Licensed under the Apache License, Version 2.0 * + * * + * Modifications: * + * - Translated from Python to C++ * + * * + ******************************************************************************/ + +#include +#include +#include +#include + +#include "cudaq/solvers/operators/molecule/fermion_compilers/bravyi_kitaev.h" + +using namespace cudaqx; + +namespace cudaq::solvers { + +template +std::set set_difference(const std::set &set1, const std::set &set2) { + std::set result; + std::set_difference(set1.begin(), set1.end(), set2.begin(), set2.end(), + std::inserter(result, result.begin())); + return result; +} + +template +std::set set_union(const std::set &set1, const std::set &set2) { + std::set result; + std::set_union(set1.begin(), set1.end(), set2.begin(), set2.end(), + std::inserter(result, result.begin())); + return result; +} + +template +std::set set_intersection(const std::set &set1, const std::set &set2) { + std::set result; + std::set_intersection(set1.begin(), set1.end(), set2.begin(), set2.end(), + std::inserter(result, result.begin())); + return result; +} + +template +std::set set_symmetric_difference(const std::set &set1, + const std::set &set2) { + std::set result; + std::set_symmetric_difference(set1.begin(), set1.end(), set2.begin(), + set2.end(), + std::inserter(result, result.begin())); + return result; +} + +/// @brief Contains the indices of qubits in the Bravyi-Kitaev basis that +/// contribute to the occupation of index-th qubit. +std::set occupation_set(std::size_t index) { + std::set indices; + index += 1; + indices.insert(index - 1); + + std::size_t parent = index & (index - 1); + index -= 1; + + while (index != parent) { + indices.insert(index - 1); + index &= index - 1; + } + + return indices; +} + +/// @brief Contains the indices of qubits in the Bravyi-Kitaev basis that +/// contribute to the parity of the index-th qubit. +std::set parity_set(std::size_t index) { + std::set indices; + + while (index > 0) { + indices.insert(index - 1); + index &= index - 1; + } + + return indices; +} + +/// @brief Contains the indices of qubits in the Bravyi-Kitaev basis that need +/// to be updated when the occupation of the index-th qubit changes. +std::set update_set(std::size_t index, std::size_t n_qubits) { + std::set indices; + + index += 1; + index += index & -index; + + while (index <= n_qubits) { + indices.insert(index - 1); + index += index & -index; + } + + return indices; +} + +std::set remainder_set(int index) { + return set_difference(parity_set(index), occupation_set(index)); +} + +std::set F_ij_set(std::size_t i, std::size_t j) { + return set_symmetric_difference(occupation_set(i), occupation_set(j)); +} + +std::set P0_ij_set(std::size_t i, std::size_t j) { + return set_symmetric_difference(parity_set(i), parity_set(j)); +} + +std::set P1_ij_set(std::size_t i, std::size_t j) { + return set_symmetric_difference(parity_set(i), remainder_set(j)); +} + +std::set P2_ij_set(std::size_t i, std::size_t j) { + return set_symmetric_difference(remainder_set(i), parity_set(j)); +} + +std::set P3_ij_set(std::size_t i, std::size_t j) { + return set_symmetric_difference(remainder_set(i), remainder_set(j)); +} + +std::set U_ij_set(std::size_t i, std::size_t j, + std::size_t n_qubits) { + return set_symmetric_difference(update_set(i, n_qubits), + update_set(j, n_qubits)); +} + +std::set alpha_set(std::size_t i, std::size_t j, + std::size_t n_qubits) { + return set_intersection(update_set(i, n_qubits), parity_set(j)); +} + +std::set U_diff_a_set(std::size_t i, std::size_t j, + std::size_t n_qubits) { + return set_difference(U_ij_set(i, j, n_qubits), alpha_set(i, j, n_qubits)); +} + +std::set P0_ij_diff_a_set(std::size_t i, std::size_t j, + std::size_t n_qubits) { + return set_symmetric_difference(P0_ij_set(i, j), alpha_set(i, j, n_qubits)); +} + +cudaq::spin_op seeley_richard_love(std::size_t i, std::size_t j, + std::complex coef, int n_qubits) { + + using double_complex = std::complex; + const double_complex imag_i = double_complex(0.0, 1.0); + + coef *= 0.25; + + cudaq::spin_op seeley_richard_love_result = 0.0 * cudaq::spin::i(0); + + // Case 0 + if (i == j) { + if (not occupation_set(i).empty()) { + cudaq::spin_op ops; + for (int index : occupation_set(i)) { + ops *= cudaq::spin::z(index); + } + seeley_richard_love_result += -coef * 2.0 * ops; + } + seeley_richard_love_result += coef * 2.0 * cudaq::spin::i(0); + } + + // Case 1 + else if (i % 2 == 0 and j % 2 == 0) { + cudaq::spin_op x_pad; + for (int index : U_diff_a_set(i, j, n_qubits)) { + x_pad *= cudaq::spin::x(index); + } + cudaq::spin_op y_pad; + for (int index : alpha_set(i, j, n_qubits)) { + y_pad *= cudaq::spin::y(index); + } + cudaq::spin_op z_pad; + for (int index : P0_ij_diff_a_set(i, j, n_qubits)) { + z_pad *= cudaq::spin::z(index); + } + cudaq::spin_op left_pad = x_pad * y_pad * z_pad; + + cudaq::spin_op op1 = left_pad * cudaq::spin::y(j) * cudaq::spin::x(i); + cudaq::spin_op op2 = left_pad * cudaq::spin::x(j) * cudaq::spin::y(i); + cudaq::spin_op op3 = left_pad * cudaq::spin::x(j) * cudaq::spin::x(i); + cudaq::spin_op op4 = left_pad * cudaq::spin::y(j) * cudaq::spin::y(i); + + if (i < j) { + seeley_richard_love_result += coef * op1; + seeley_richard_love_result += -coef * op2; + seeley_richard_love_result += -imag_i * coef * op3; + seeley_richard_love_result += -imag_i * coef * op4; + } else { + seeley_richard_love_result += -imag_i * coef * op1; + seeley_richard_love_result += imag_i * coef * op2; + seeley_richard_love_result += -coef * op3; + seeley_richard_love_result += -coef * op4; + } + } + + // Case 2 + else if (i % 2 == 1 and j % 2 == 0 and not parity_set(j).contains(i)) { + cudaq::spin_op x_pad; + for (int index : U_diff_a_set(i, j, n_qubits)) { + x_pad *= cudaq::spin::x(index); + } + cudaq::spin_op y_pad; + for (int index : alpha_set(i, j, n_qubits)) { + y_pad *= cudaq::spin::y(index); + } + cudaq::spin_op left_pad = x_pad * y_pad; + + cudaq::spin_op right_pad_1; + auto P0_minus_alpha = + set_difference(P0_ij_set(i, j), alpha_set(i, j, n_qubits)); + for (int index : P0_minus_alpha) { + right_pad_1 *= cudaq::spin::z(index); + } + + cudaq::spin_op right_pad_2; + auto P2_minus_alpha = + set_difference(P2_ij_set(i, j), alpha_set(i, j, n_qubits)); + for (int index : P2_minus_alpha) { + right_pad_2 *= cudaq::spin::z(index); + } + + double_complex c0, c1, c2, c3; + if (i < j) { + c0 = coef; + c1 = -imag_i * coef; + c2 = -coef; + c3 = -imag_i * coef; + } else { + c0 = -imag_i * coef; + c1 = -coef; + c2 = imag_i * coef; + c3 = -coef; + } + + seeley_richard_love_result += + c0 * left_pad * cudaq::spin::y(j) * cudaq::spin::x(i) * right_pad_1; + seeley_richard_love_result += + c1 * left_pad * cudaq::spin::x(j) * cudaq::spin::x(i) * right_pad_1; + seeley_richard_love_result += + c2 * left_pad * cudaq::spin::x(j) * cudaq::spin::y(i) * right_pad_2; + seeley_richard_love_result += + c3 * left_pad * cudaq::spin::y(j) * cudaq::spin::y(i) * right_pad_2; + } + + // Case 3 + else if (i % 2 == 1 and j % 2 == 0 and parity_set(j).contains(i)) { + cudaq::spin_op left_pad; + for (auto index : U_ij_set(i, j, n_qubits)) { + left_pad *= cudaq::spin::x(index); + } + + cudaq::spin_op right_pad_1; + for (auto index : set_difference(P0_ij_set(i, j), {i})) { + right_pad_1 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_2; + for (auto index : set_difference(P2_ij_set(i, j), {i})) { + right_pad_2 *= cudaq::spin::z(index); + } + + seeley_richard_love_result += + coef * left_pad * cudaq::spin::y(j) * cudaq::spin::y(i) * right_pad_1; + seeley_richard_love_result += -imag_i * coef * left_pad * + cudaq::spin::x(j) * cudaq::spin::y(i) * + right_pad_1; + seeley_richard_love_result += + coef * left_pad * cudaq::spin::x(j) * cudaq::spin::x(i) * right_pad_2; + seeley_richard_love_result += imag_i * coef * left_pad * cudaq::spin::y(j) * + cudaq::spin::x(i) * right_pad_2; + } + + // Case 4 + else if (i % 2 == 0 and j % 2 == 1 and not parity_set(j).contains(i) and + not update_set(i, n_qubits).contains(i)) { + cudaq::spin_op x_pad; + for (auto index : U_diff_a_set(i, j, n_qubits)) { + x_pad *= cudaq::spin::x(index); + } + cudaq::spin_op y_pad; + for (auto index : alpha_set(i, j, n_qubits)) { + y_pad *= cudaq::spin::y(index); + } + auto left_pad = x_pad * y_pad; + + cudaq::spin_op right_pad_1; + for (auto index : + set_difference(P0_ij_set(i, j), alpha_set(i, j, n_qubits))) { + right_pad_1 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_2; + for (auto index : + set_difference(P1_ij_set(i, j), alpha_set(i, j, n_qubits))) { + right_pad_2 *= cudaq::spin::z(index); + } + + double_complex c0, c1, c2, c3; + if (i < j) { + c0 = -coef; + c1 = -imag_i * coef; + c2 = coef; + c3 = -imag_i * coef; + } else { + c0 = imag_i * coef; + c1 = -coef; + c2 = -imag_i * coef; + c3 = -coef; + } + seeley_richard_love_result += + c0 * left_pad * cudaq::spin::x(j) * cudaq::spin::y(i) * right_pad_1; + seeley_richard_love_result += + c1 * left_pad * cudaq::spin::x(j) * cudaq::spin::x(i) * right_pad_1; + seeley_richard_love_result += + c2 * left_pad * cudaq::spin::y(j) * cudaq::spin::x(i) * right_pad_2; + seeley_richard_love_result += + c3 * left_pad * cudaq::spin::y(j) * cudaq::spin::y(i) * right_pad_2; + } + + // Case 5 + else if (i % 2 == 0 and j % 2 == 1 and not parity_set(j).contains(i) and + update_set(i, n_qubits).contains(j)) { + cudaq::spin_op left_pad_1; + auto x_range_1 = set_difference(U_ij_set(i, j, n_qubits), {j}); + for (auto index : x_range_1) { + left_pad_1 *= cudaq::spin::x(index); + } + cudaq::spin_op left_pad_2; + auto x_range_2 = set_difference(x_range_1, alpha_set(i, j, n_qubits)); + for (auto index : x_range_2) { + left_pad_2 *= cudaq::spin::x(index); + } + + cudaq::spin_op y_pad; + for (auto index : alpha_set(i, j, n_qubits)) { + y_pad *= cudaq::spin::y(index); + } + cudaq::spin_op z_pad; + for (auto index : + set_difference(P0_ij_set(i, j), alpha_set(i, j, n_qubits))) { + z_pad *= cudaq::spin::z(index); + } + auto right_pad_1 = y_pad * z_pad; + + cudaq::spin_op right_pad_2; + for (auto index : set_union(P1_ij_set(i, j), {j})) { + right_pad_2 *= cudaq::spin::z(index); + } + } + + // Case 6 + else if (i % 2 == 0 and j % 2 == 1 and parity_set(j).contains(i) and + update_set(i, n_qubits).contains(j)) { + cudaq::spin_op left_pad; + for (auto index : set_difference(U_ij_set(i, j, n_qubits), {j})) { + left_pad *= cudaq::spin::x(index); + } + cudaq::spin_op right_pad; + for (auto index : set_union(P1_ij_set(i, j), {j})) { + right_pad *= cudaq::spin::z(index); + } + + seeley_richard_love_result += coef * left_pad * cudaq::spin::x(i); + seeley_richard_love_result += -imag_i * coef * left_pad * cudaq::spin::y(i); + seeley_richard_love_result += + imag_i * coef * left_pad * cudaq::spin::y(i) * right_pad; + seeley_richard_love_result += + -coef * left_pad * cudaq::spin::x(i) * right_pad; + } + + // Case 7 + else if (i % 2 == 1 and j % 2 == 1 and not parity_set(j).contains(i) and + not update_set(i, n_qubits).contains(j)) { + cudaq::spin_op x_pad; + for (auto index : U_diff_a_set(i, j, n_qubits)) { + x_pad *= cudaq::spin::x(index); + } + cudaq::spin_op y_pad; + for (auto index : alpha_set(i, j, n_qubits)) { + y_pad *= cudaq::spin::y(index); + } + auto left_pad = x_pad * y_pad; + + cudaq::spin_op right_pad_1; + for (auto index : + set_difference(P0_ij_set(i, j), alpha_set(i, j, n_qubits))) { + right_pad_1 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_2; + for (auto index : + set_difference(P1_ij_set(i, j), alpha_set(i, j, n_qubits))) { + right_pad_2 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_3; + for (auto index : + set_difference(P2_ij_set(i, j), alpha_set(i, j, n_qubits))) { + right_pad_3 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_4; + for (auto index : + set_difference(P3_ij_set(i, j), alpha_set(i, j, n_qubits))) { + right_pad_4 *= cudaq::spin::z(index); + } + + double_complex c0, c1, c2, c3; + if (i < j) { + c0 = -imag_i * coef; + c1 = coef; + c2 = -coef; + c3 = -imag_i * coef; + } else { + c0 = -coef; + c1 = -imag_i * coef; + c2 = imag_i * coef; + c3 = -coef; + } + + seeley_richard_love_result += + c0 * left_pad * cudaq::spin::x(j) * cudaq::spin::x(i) * right_pad_1; + seeley_richard_love_result += + c1 * left_pad * cudaq::spin::y(j) * cudaq::spin::x(i) * right_pad_2; + seeley_richard_love_result += + c2 * left_pad * cudaq::spin::x(j) * cudaq::spin::y(i) * right_pad_3; + seeley_richard_love_result += + c3 * left_pad * cudaq::spin::y(j) * cudaq::spin::y(i) * right_pad_4; + } + + // Case 8 + else if (i % 2 == 1 and j % 2 == 1 and parity_set(j).contains(i) and + not update_set(i, n_qubits).contains(j)) { + cudaq::spin_op left_pad; + for (auto index : U_ij_set(i, j, n_qubits)) { + left_pad *= cudaq::spin::x(index); + } + + cudaq::spin_op right_pad_1; + for (auto index : set_difference(P0_ij_set(i, j), {i})) { + right_pad_1 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_2; + for (auto index : set_difference(P1_ij_set(i, j), {i})) { + right_pad_2 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_3; + for (auto index : set_difference(P2_ij_set(i, j), {i})) { + right_pad_3 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_4; + for (auto index : set_difference(P3_ij_set(i, j), {i})) { + right_pad_4 *= cudaq::spin::z(index); + } + + seeley_richard_love_result += -imag_i * coef * left_pad * + cudaq::spin::x(j) * cudaq::spin::y(i) * + right_pad_1; + seeley_richard_love_result += + coef * left_pad * cudaq::spin::y(j) * cudaq::spin::y(i) * right_pad_2; + seeley_richard_love_result += + coef * left_pad * cudaq::spin::x(j) * cudaq::spin::x(i) * right_pad_3; + seeley_richard_love_result += imag_i * coef * left_pad * cudaq::spin::y(j) * + cudaq::spin::x(i) * right_pad_4; + } + + // Case 9 + else if (i % 2 == 1 and j % 2 == 1 and not parity_set(j).contains(i) and + update_set(i, n_qubits).contains(j)) { + cudaq::spin_op left_pad_3; + auto x_range_1 = set_difference(U_ij_set(i, j, n_qubits), {j}); + for (auto index : x_range_1) { + left_pad_3 *= cudaq::spin::x(index); + } + + cudaq::spin_op left_pad_1; + auto x_range_2 = set_difference(x_range_1, alpha_set(i, j, n_qubits)); + for (auto index : x_range_2) { + left_pad_1 *= cudaq::spin::x(index); + } + + cudaq::spin_op left_pad_2; + auto x_range_3 = + set_difference(set_union(x_range_1, {i}), alpha_set(i, j, n_qubits)); + for (auto index : x_range_3) { + left_pad_2 *= cudaq::spin::x(index); + } + + auto z_range_1 = set_difference(P2_ij_set(i, j), alpha_set(i, j, n_qubits)); + cudaq::spin_op z_pad_1; + for (auto index : z_range_1) { + z_pad_1 *= cudaq::spin::z(index); + } + cudaq::spin_op y_pad; + for (auto index : alpha_set(i, j, n_qubits)) { + y_pad *= cudaq::spin::y(index); + } + auto right_pad_1 = z_pad_1 * y_pad; + + auto z_range_2 = set_difference(P0_ij_set(i, j), alpha_set(i, j, n_qubits)); + cudaq::spin_op z_pad_2; + for (auto index : z_range_2) { + z_pad_2 *= cudaq::spin::z(index); + } + auto right_pad_2 = z_pad_2 * y_pad; + + auto z_range_3 = set_union(P1_ij_set(i, j), {j}); + cudaq::spin_op right_pad_3; + for (auto index : z_range_3) { + right_pad_3 *= cudaq::spin::z(index); + } + + auto z_range_4 = set_union(P3_ij_set(i, j), {j}); + cudaq::spin_op right_pad_4; + for (auto index : z_range_4) { + right_pad_4 *= cudaq::spin::z(index); + } + + seeley_richard_love_result += + -coef * left_pad_1 * cudaq::spin::y(i) * right_pad_1; + seeley_richard_love_result += -imag_i * coef * left_pad_2 * right_pad_2; + seeley_richard_love_result += + -coef * left_pad_3 * cudaq::spin::x(i) * right_pad_3; + seeley_richard_love_result += + imag_i * coef * left_pad_3 * cudaq::spin::y(i) * right_pad_4; + } + + // Case 10 + else if (i % 2 == 1 and j % 2 == 1 and parity_set(j).contains(i) and + update_set(i, n_qubits).contains(j)) { + cudaq::spin_op left_pad; + for (auto index : set_difference(U_ij_set(i, j, n_qubits), {j})) { + left_pad *= cudaq::spin::x(index); + } + cudaq::spin_op right_pad_1; + for (auto index : set_difference(P0_ij_set(i, j), {i})) { + right_pad_1 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_2; + for (auto index : set_difference(P2_ij_set(i, j), {i})) { + right_pad_2 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_3; + for (auto index : P1_ij_set(i, j)) { + right_pad_3 *= cudaq::spin::z(index); + } + cudaq::spin_op right_pad_4; + for (auto index : P3_ij_set(i, j)) { + right_pad_4 *= cudaq::spin::z(index); + } + + seeley_richard_love_result += + -imag_i * coef * left_pad * cudaq::spin::y(i) * right_pad_1; + seeley_richard_love_result += + coef * left_pad * cudaq::spin::x(i) * right_pad_2; + seeley_richard_love_result += + -coef * left_pad * cudaq::spin::z(j) * cudaq::spin::x(i) * right_pad_3; + seeley_richard_love_result += imag_i * coef * left_pad * cudaq::spin::z(j) * + cudaq::spin::y(i) * right_pad_4; + } + + return seeley_richard_love_result; +} + +cudaq::spin_op hermitian_one_body_product(std::size_t a, std::size_t b, + std::size_t c, std::size_t d, + std::complex coef, + std::size_t nqubits) { + + cudaq::spin_op c_dag_c_ac = seeley_richard_love(a, c, coef, nqubits); + cudaq::spin_op c_dag_c_bd = seeley_richard_love(b, d, 1, nqubits); + c_dag_c_ac *= c_dag_c_bd; + + cudaq::spin_op hermitian_sum = c_dag_c_ac; + + cudaq::spin_op c_dag_c_ca = + seeley_richard_love(c, a, std::conj(coef), nqubits); + cudaq::spin_op c_dag_c_db = seeley_richard_love(d, b, 1, nqubits); + c_dag_c_ca *= c_dag_c_db; + + hermitian_sum += c_dag_c_ca; + return hermitian_sum; +} + +std::complex two_body_coef(const cudaqx::tensor<> &hpqrs, std::size_t p, + std::size_t q, std::size_t r, + std::size_t s) { + return hpqrs.at({p, q, r, s}) - hpqrs.at({q, p, r, s}) - + hpqrs.at({p, q, s, r}) + hpqrs.at({q, p, s, r}); +} + +cudaq::spin_op bravyi_kitaev::generate(const double constant, + const tensor<> &hpq, + const tensor<> &hpqrs, + const heterogeneous_map &options) { + assert(hpq.rank() == 2 && "hpq must be a rank-2 tensor"); + assert(hpqrs.rank() == 4 && "hpqrs must be a rank-4 tensor"); + auto nqubits = hpq.shape()[0]; + cudaq::spin_op spin_hamiltonian = 0.0 * cudaq::spin::i(0); + + double tolerance = + options.get(std::vector{"tolerance", "tol"}, 1e-15); + + double constant_term = constant; + for (std::size_t p = 0; p < nqubits; p++) { + if (std::abs(hpq.at({p, p})) > 0.0) { + spin_hamiltonian += seeley_richard_love(p, p, hpq.at({p, p}), nqubits); + } + for (std::size_t q = 0; q < p; q++) { + if (std::abs(hpq.at({p, q})) > 0.0) { + spin_hamiltonian += seeley_richard_love(p, q, hpq.at({p, q}), nqubits); + spin_hamiltonian += + seeley_richard_love(q, p, std::conj(hpq.at({p, q})), nqubits); + } + auto coef = 0.25 * two_body_coef(hpqrs, p, q, q, p); + if (std::abs(coef) > 0.0) { + if (not occupation_set(p).empty()) { + cudaq::spin_op zs; + for (auto index : occupation_set(p)) { + zs *= cudaq::spin::z(index); + } + spin_hamiltonian += -coef * zs; + } + if (not occupation_set(q).empty()) { + cudaq::spin_op zs2; + for (auto index : occupation_set(q)) { + zs2 *= cudaq::spin::z(index); + } + spin_hamiltonian += -coef * zs2; + } + if (not F_ij_set(p, q).empty()) { + cudaq::spin_op zs3; + for (auto index : F_ij_set(p, q)) { + zs3 *= cudaq::spin::z(index); + } + spin_hamiltonian += coef * zs3; + } + constant_term += std::real(coef); + } + } + } + + for (std::size_t p = 0; p < nqubits; p++) { + for (std::size_t q = 0; q < nqubits; q++) { + for (std::size_t r = 0; r < q; r++) { + if ((p != q) and (p != r)) { + auto coef = two_body_coef(hpqrs, p, q, r, p); + if (std::abs(coef) > 0.0) { + cudaq::spin_op excitation = + seeley_richard_love(q, r, coef, nqubits); + cudaq::spin_op number = seeley_richard_love(p, p, 1.0, nqubits); + spin_hamiltonian += number * excitation; + } + } + } + } + } + + for (std::size_t p = 0; p < nqubits; p++) { + for (std::size_t q = 0; q < p; q++) { + for (std::size_t r = 0; r < q; r++) { + for (std::size_t s = 0; s < r; s++) { + auto coef_pqrs = -two_body_coef(hpqrs, p, q, r, s); + if (std::abs(coef_pqrs) > 0.0) { + spin_hamiltonian += + hermitian_one_body_product(p, q, r, s, coef_pqrs, nqubits); + } + auto coef_prqs = -two_body_coef(hpqrs, p, r, q, s); + if (std::abs(coef_prqs) > 0.0) { + spin_hamiltonian += + hermitian_one_body_product(p, r, q, s, coef_prqs, nqubits); + } + auto coef_psqr = -two_body_coef(hpqrs, p, s, q, r); + if (std::abs(coef_psqr) > 0.0) { + spin_hamiltonian += + hermitian_one_body_product(p, s, q, r, coef_psqr, nqubits); + } + } + } + } + } + + spin_hamiltonian += constant_term * cudaq::spin::i(0); + return spin_hamiltonian; +} +} // namespace cudaq::solvers diff --git a/libs/solvers/python/bindings/solvers/py_solvers.cpp b/libs/solvers/python/bindings/solvers/py_solvers.cpp index 16aa7b3..fbc88aa 100644 --- a/libs/solvers/python/bindings/solvers/py_solvers.cpp +++ b/libs/solvers/python/bindings/solvers/py_solvers.cpp @@ -353,6 +353,158 @@ RuntimeError further manipulated using CUDA Quantum operations. )#"); + mod.def( + "bravyi_kitaev", + [](py::buffer hpq, py::buffer hpqrs, double core_energy = 0.0, + py::kwargs options) { + auto hpqInfo = hpq.request(); + auto hpqrsInfo = hpqrs.request(); + auto *hpqData = reinterpret_cast *>(hpqInfo.ptr); + auto *hpqrsData = + reinterpret_cast *>(hpqrsInfo.ptr); + + cudaqx::tensor hpqT, hpqrsT; + hpqT.borrow(hpqData, {hpqInfo.shape.begin(), hpqInfo.shape.end()}); + hpqrsT.borrow(hpqrsData, + {hpqrsInfo.shape.begin(), hpqrsInfo.shape.end()}); + + return fermion_compiler::get("bravyi_kitaev") + ->generate(core_energy, hpqT, hpqrsT, hetMapFromKwargs(options)); + }, + py::arg("hpq"), py::arg("hpqrs"), py::arg("core_energy") = 0.0, + R"#( +Perform the Bravyi-Kitaev transformation on fermionic operators. + +This function applies the Bravyi-Kitaev transformation to convert fermionic operators +(represented by one- and two-body integrals) into qubit operators. + +Parameters: +----------- +hpq : numpy.ndarray + A 2D complex numpy array representing the one-body integrals. + Shape should be (N, N) where N is the number of spin molecular orbitals. +hpqrs : numpy.ndarray + A 4D complex numpy array representing the two-body integrals. + Shape should be (N, N, N, N) where N is the number of spin molecular orbitals. +core_energy : float, optional + The core energy of the system when using active space Hamiltonian, nuclear energy otherwise. Default is 0.0. +tolerance : float, optional + The threshold value for ignoring small coefficients. + Can also be specified using 'tol'. + Coefficients with absolute values smaller than this tolerance are considered as zero. + Default is 1e-15. + +Returns: +-------- +cudaq.SpinOperator + A qubit operator (spin operator) resulting from the Bravyi-Kitaev transformation. + +Raises: +------- +ValueError + If the input arrays have incorrect shapes or types. +RuntimeError + If the Bravyi-Kitaev transformation fails for any reason. + +Examples: +--------- +>>> import numpy as np +>>> h1 = np.array([[0, 1], [1, 0]], dtype=np.complex128) +>>> h2 = np.zeros((2, 2, 2, 2), dtype=np.complex128) +>>> h2[0, 1, 1, 0] = h2[1, 0, 0, 1] = 0.5 +>>> qubit_op = bravyi_kitaev(h1, h2, core_energy=0.1, tolerance=1e-14) + +Notes: +------ +- The input arrays `hpq` and `hpqrs` must be contiguous and in row-major order. +- This function uses the "bravyi_kitaev" fermion compiler internally to perform + the transformation. +- The resulting qubit operator can be used directly in quantum algorithms or + further manipulated using CUDA Quantum operations. +)#"); + + mod.def( + "bravyi_kitaev", + [](py::buffer buffer, double core_energy = 0.0, py::kwargs options) { + auto info = buffer.request(); + auto *data = reinterpret_cast *>(info.ptr); + std::size_t size = 1; + for (auto &s : info.shape) + size *= s; + std::vector> vec(data, data + size); + if (info.shape.size() == 2) { + std::size_t dim = info.shape[0]; + cudaqx::tensor hpq, hpqrs({dim, dim, dim, dim}); + hpq.borrow(data, {info.shape.begin(), info.shape.end()}); + return fermion_compiler::get("bravyi_kitaev") + ->generate(core_energy, hpq, hpqrs, hetMapFromKwargs(options)); + } + + std::size_t dim = info.shape[0]; + cudaqx::tensor hpq({dim, dim}), hpqrs; + hpqrs.borrow(data, {info.shape.begin(), info.shape.end()}); + return fermion_compiler::get("bravyi_kitaev") + ->generate(core_energy, hpq, hpqrs, hetMapFromKwargs(options)); + }, + py::arg("hpq"), py::arg("core_energy") = 0.0, + R"#( +Perform the Bravyi-Kitaev transformation on fermionic operators. + +This function applies the Bravyi-Kitaev transformation to convert fermionic operators +(represented by either one-body or two-body integrals) into qubit operators. + +Parameters: +----------- +hpq : numpy.ndarray + A complex numpy array representing either: + - One-body integrals: 2D array with shape (N, N) + - Two-body integrals: 4D array with shape (N, N, N, N) + where N is the number of orbitals. +core_energy : float, optional + The core energy of the system. Default is 0.0. +tolerance : float, optional + The threshold value for ignoring small coefficients. + Can also be specified using 'tol'. + Coefficients with absolute values smaller than this tolerance are considered as zero. + Default is 1e-15. + +Returns: +-------- +cudaq.SpinOperator + A qubit operator (spin operator) resulting from the Bravyi-Kitaev transformation. + +Raises: +------- +ValueError + If the input array has an incorrect shape or type. +RuntimeError + If the Bravyi-Kitaev transformation fails for any reason. + +Examples: +--------- +>>> import numpy as np +>>> # One-body integrals +>>> h1 = np.array([[0, 1], [1, 0]], dtype=np.complex128) +>>> qubit_op1 = bravyi_kitaev(h1, core_energy=0.1, tolerance=1e-14) + +>>> # Two-body integrals +>>> h2 = np.zeros((2, 2, 2, 2), dtype=np.complex128) +>>> h2[0, 1, 1, 0] = h2[1, 0, 0, 1] = 0.5 +>>> qubit_op2 = bravyi_kitaev(h2) + +Notes: +------ +- The input array must be contiguous and in row-major order. +- This function automatically detects whether the input represents one-body or + two-body integrals based on its shape. +- For one-body integrals input, a zero-initialized two-body tensor is used internally. +- For two-body integrals input, a zero-initialized one-body tensor is used internally. +- This function uses the "bravyi_kitaev" fermion compiler internally to perform + the transformation. +- The resulting qubit operator can be used directly in quantum algorithms or + further manipulated using CUDA Quantum operations. +)#"); + py::class_(mod, "MolecularHamiltonian") .def_readonly("energies", &molecular_hamiltonian::energies, R"#( diff --git a/libs/solvers/python/tests/test_molecule.py b/libs/solvers/python/tests/test_molecule.py index df837c0..041b773 100644 --- a/libs/solvers/python/tests/test_molecule.py +++ b/libs/solvers/python/tests/test_molecule.py @@ -75,6 +75,28 @@ def test_jordan_wigner(): assert np.isclose(np.min(e), -1.13717, rtol=1e-4) +def test_bravyi_kitaev(): + geometry = [('H', (0., 0., 0.)), ('H', (0., 0., .7474))] + molecule = solvers.create_molecule(geometry, + 'sto-3g', + 0, + 0, + verbose=True, + casci=True) + + op = solvers.bravyi_kitaev(molecule.hpq, + molecule.hpqrs, + core_energy=molecule.energies['nuclear_energy'], + tol=1e-15) + spin_ham_matrix = molecule.hamiltonian.to_matrix() + e, c = np.linalg.eig(spin_ham_matrix) + assert np.isclose(np.min(e), -1.13717, rtol=1e-4) + + spin_ham_matrix = op.to_matrix() + e, c = np.linalg.eig(spin_ham_matrix) + assert np.isclose(np.min(e), -1.13717, rtol=1e-4) + + def test_active_space(): geometry = [('N', (0.0, 0.0, 0.5600)), ('N', (0.0, 0.0, -0.5600))] @@ -135,6 +157,31 @@ def test_jordan_wigner_as(): assert np.isclose(np.min(e), -107.542198, rtol=1e-4) +def test_bravyi_kitaev_as(): + geometry = [('N', (0.0, 0.0, 0.5600)), ('N', (0.0, 0.0, -0.5600))] + molecule = solvers.create_molecule(geometry, + 'sto-3g', + 0, + 0, + nele_cas=4, + norb_cas=4, + ccsd=True, + casci=True, + verbose=True) + + op = solvers.bravyi_kitaev(molecule.hpq, + molecule.hpqrs, + core_energy=molecule.energies['core_energy'], + tol=1e-15) + spin_ham_matrix = molecule.hamiltonian.to_matrix() + e, c = np.linalg.eig(spin_ham_matrix) + assert np.isclose(np.min(e), -107.542198, rtol=1e-4) + + spin_ham_matrix = op.to_matrix() + e, c = np.linalg.eig(spin_ham_matrix) + assert np.isclose(np.min(e), -107.542198, rtol=1e-4) + + def test_as_with_natorb(): geometry = [('N', (0.0, 0.0, 0.5600)), ('N', (0.0, 0.0, -0.5600))] diff --git a/libs/solvers/unittests/CMakeLists.txt b/libs/solvers/unittests/CMakeLists.txt index 2982583..4716311 100644 --- a/libs/solvers/unittests/CMakeLists.txt +++ b/libs/solvers/unittests/CMakeLists.txt @@ -44,6 +44,11 @@ target_link_libraries(test_molecule PRIVATE GTest::gtest_main cudaq-solvers) add_dependencies(CUDAQXSolversUnitTests test_molecule) gtest_discover_tests(test_molecule) +add_executable(test_bravyi_kitaev test_bravyi_kitaev.cpp) +target_link_libraries(test_bravyi_kitaev PRIVATE GTest::gtest_main cudaq-solvers) +add_dependencies(CUDAQXSolversUnitTests test_bravyi_kitaev) +gtest_discover_tests(test_bravyi_kitaev) + add_executable(test_optimizers test_optimizers.cpp) target_link_libraries(test_optimizers PRIVATE GTest::gtest_main cudaq-solvers) add_dependencies(CUDAQXSolversUnitTests test_optimizers) @@ -66,4 +71,4 @@ gtest_discover_tests(test_uccsd) add_executable(test_qaoa test_qaoa.cpp) target_link_libraries(test_qaoa PRIVATE GTest::gtest_main cudaq-solvers test-kernels) -gtest_discover_tests(test_qaoa) \ No newline at end of file +gtest_discover_tests(test_qaoa) diff --git a/libs/solvers/unittests/support/h2_pyscf_hf.py b/libs/solvers/unittests/support/h2_pyscf_hf.py new file mode 100644 index 0000000..747302c --- /dev/null +++ b/libs/solvers/unittests/support/h2_pyscf_hf.py @@ -0,0 +1,73 @@ +""" +H2 Hamiltonian Generator +======================= + +This script generates ground truth data for the test_bravyi_kitaev.cpp unit tests by creating +and transforming a molecular Hamiltonian for the H2 molecule. The process involves: + +1. Creating a molecular Hamiltonian using PySCF (via OpenFermion) +2. Extracting one- and two-body integrals to be used as coefficients +3. Converting to a spin Hamiltonian using the Bravyi-Kitaev transformation + +The output includes: +- System information (number of orbitals, electrons) +- Nuclear repulsion energy +- Hartree-Fock energy +- Molecular Hamiltonian terms +- Qubit Hamiltonian terms after Bravyi-Kitaev transformation + +Dependencies: + - numpy + - openfermion + - openfermionpyscf + +Configuration: + - Molecule: H2 + - Geometry: H atoms at [0,0,0] and [0,0,0.7474] + - Basis set: STO-3G + - Multiplicity: 1 (singlet) + - Charge: 0 (neutral) +""" + +import numpy as np +from openfermion import * +from openfermionpyscf import run_pyscf + +geometry = [['H', [0, 0, 0]], ['H', [0, 0, 0.7474]]] +basis = 'sto-3g' +multiplicity = 1 +charge = 0 +molecule = MolecularData(geometry, basis, multiplicity, charge) +molecule = run_pyscf(molecule, run_scf=True) + +# Get the molecular Hamiltonian +molecular_hamiltonian = molecule.get_molecular_hamiltonian() + +# Convert to fermion operator +fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian) + +# Convert to qubit Hamiltonian using Bravyi-Kitaev transformation +qubit_hamiltonian = bravyi_kitaev(fermion_hamiltonian) + + +def print_hamiltonian_info(): + print("System Information:") + print(f"Number of orbitals: {molecule.n_orbitals}") + print(f"Number of electrons: {molecule.n_electrons}") + print(f"Nuclear repulsion energy: {molecule.nuclear_repulsion:.8f}") + print(f"HF energy: {molecule.hf_energy:.8f}") + + print("\nMolecular Hamiltonian terms:") + print(molecular_hamiltonian) + + print("\nNumber of qubits required:") + print(count_qubits(qubit_hamiltonian)) + + print("\nQubit Hamiltonian terms:") + for term, coefficient in qubit_hamiltonian.terms.items(): + if abs(coefficient) > 1e-8: # Filter out near-zero terms + print(f"{coefficient:.8f} [{' '.join(str(x) for x in term)}]") + + +# Print the Hamiltonians and system information +print_hamiltonian_info() diff --git a/libs/solvers/unittests/test_bravyi_kitaev.cpp b/libs/solvers/unittests/test_bravyi_kitaev.cpp new file mode 100644 index 0000000..8f95e2c --- /dev/null +++ b/libs/solvers/unittests/test_bravyi_kitaev.cpp @@ -0,0 +1,243 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2024 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ +#include +#include +#include +#include +#include + +#include + +#include "cudaq/solvers/operators/molecule/fermion_compilers/bravyi_kitaev.h" + +// One- and Two-body integrals were copied from test_molecule.cpp. +// They were further validated using the script ./support/h2_pyscf_hf.py. +// +TEST(BravyiKitaev, testH2Hamiltonian) { + using double_complex = std::complex; + using namespace cudaq::spin; + + cudaqx::tensor<> hpq({4, 4}); + cudaqx::tensor<> hpqrs({4, 4, 4, 4}); + + double h_constant = 0.7080240981000804; + hpq.at({0, 0}) = -1.2488; + hpq.at({1, 1}) = -1.2488; + hpq.at({2, 2}) = -.47967; + hpq.at({3, 3}) = -.47967; + hpqrs.at({0, 0, 0, 0}) = 0.3366719725032414; + hpqrs.at({0, 0, 2, 2}) = 0.0908126657382825; + hpqrs.at({0, 1, 1, 0}) = 0.3366719725032414; + hpqrs.at({0, 1, 3, 2}) = 0.0908126657382825; + hpqrs.at({0, 2, 0, 2}) = 0.09081266573828267; + hpqrs.at({0, 2, 2, 0}) = 0.33121364716348484; + hpqrs.at({0, 3, 1, 2}) = 0.09081266573828267; + hpqrs.at({0, 3, 3, 0}) = 0.33121364716348484; + hpqrs.at({1, 0, 0, 1}) = 0.3366719725032414; + hpqrs.at({1, 0, 2, 3}) = 0.0908126657382825; + hpqrs.at({1, 1, 1, 1}) = 0.3366719725032414; + hpqrs.at({1, 1, 3, 3}) = 0.0908126657382825; + hpqrs.at({1, 2, 0, 3}) = 0.09081266573828267; + hpqrs.at({1, 2, 2, 1}) = 0.33121364716348484; + hpqrs.at({1, 3, 1, 3}) = 0.09081266573828267; + hpqrs.at({1, 3, 3, 1}) = 0.33121364716348484; + hpqrs.at({2, 0, 0, 2}) = 0.3312136471634851; + hpqrs.at({2, 0, 2, 0}) = 0.09081266573828246; + hpqrs.at({2, 1, 1, 2}) = 0.3312136471634851; + hpqrs.at({2, 1, 3, 0}) = 0.09081266573828246; + hpqrs.at({2, 2, 0, 0}) = 0.09081266573828264; + hpqrs.at({2, 2, 2, 2}) = 0.34814578499360427; + hpqrs.at({2, 3, 1, 0}) = 0.09081266573828264; + hpqrs.at({2, 3, 3, 2}) = 0.34814578499360427; + hpqrs.at({3, 0, 0, 3}) = 0.3312136471634851; + hpqrs.at({3, 0, 2, 1}) = 0.09081266573828246; + hpqrs.at({3, 1, 1, 3}) = 0.3312136471634851; + hpqrs.at({3, 1, 3, 1}) = 0.09081266573828246; + hpqrs.at({3, 2, 0, 1}) = 0.09081266573828264; + hpqrs.at({3, 2, 2, 3}) = 0.34814578499360427; + hpqrs.at({3, 3, 1, 1}) = 0.09081266573828264; + hpqrs.at({3, 3, 3, 3}) = 0.34814578499360427; + + cudaq::solvers::bravyi_kitaev transform{}; + cudaq::spin_op result = transform.generate(h_constant, hpq, hpqrs, {}); + cudaq::spin_op gold = + -0.1064770114930045 * i(0) + 0.04540633286914125 * x(0) * z(1) * x(2) + + 0.04540633286914125 * x(0) * z(1) * x(2) * z(3) + + 0.04540633286914125 * y(0) * z(1) * y(2) + + 0.04540633286914125 * y(0) * z(1) * y(2) * z(3) + + 0.17028010135220506 * z(0) + 0.1702801013522051 * z(0) * z(1) + + 0.16560682358174256 * z(0) * z(1) * z(2) + + 0.16560682358174256 * z(0) * z(1) * z(2) * z(3) + + 0.12020049071260128 * z(0) * z(2) + + 0.12020049071260128 * z(0) * z(2) * z(3) + 0.1683359862516207 * z(1) - + 0.22004130022421792 * z(1) * z(2) * z(3) + + 0.17407289249680227 * z(1) * z(3) - 0.22004130022421792 * z(2); + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase0) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(2, 2, 4.0, 20); + + cudaq::spin_op gold = double_complex(-2.0, 0.0) * i(0) * i(1) * z(2) + + double_complex(2.0, 0.0) * i(0) * i(1) * i(2); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase1) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(2, 6, 4.0, 20); + cudaq::spin_op gold = double_complex(1.0, 0.0) * i(0) * z(1) * x(2) * y(3) * + i(4) * z(5) * y(6) + + double_complex(-1.0, 0.0) * i(0) * z(1) * y(2) * y(3) * + i(4) * z(5) * x(6) + + double_complex(0.0, -1.0) * i(0) * z(1) * x(2) * y(3) * + i(4) * z(5) * x(6) + + double_complex(0.0, -1.0) * i(0) * z(1) * y(2) * y(3) * + i(4) * z(5) * y(6); + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase2) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(5, 2, 4.0, 20); + cudaq::spin_op gold = + double_complex(-1.0, 0.0) * z(1) * y(2) * y(3) * z(4) * x(5) + + double_complex(0.0, 1.0) * z(1) * x(2) * y(3) * z(4) * x(5) + + double_complex(1.0, 0.0) * z(1) * x(2) * y(3) * y(5) + + double_complex(0.0, 1.0) * z(1) * y(2) * y(3) * y(5); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase3) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(1, 2, 4.0, 20); + cudaq::spin_op gold = double_complex(1.0, 0.0) * z(0) * y(1) * y(2) + + double_complex(0.0, -1.0) * z(0) * y(1) * x(2) + + double_complex(1.0, 0.0) * i(0) * x(1) * x(2) + + double_complex(0.0, 1.0) * i(0) * x(1) * y(2); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase4) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(0, 5, 4.0, 20); + cudaq::spin_op gold = + double_complex(-1.0, 0.0) * y(0) * x(1) * i(2) * y(3) * z(4) * x(5) + + double_complex(0.0, -1.0) * x(0) * x(1) * i(2) * y(3) * z(4) * x(5) + + double_complex(1.0, 0.0) * x(0) * x(1) * i(2) * y(3) * i(4) * y(5) + + double_complex(0.0, -1.0) * y(0) * x(1) * i(2) * y(3) * i(4) * y(5); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase6) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(18, 19, 4.0, 20); + cudaq::spin_op gold = double_complex(1.0, 0.0) * x(18) * i(19) + + double_complex(0.0, -1.0) * y(18) * i(19) + + double_complex(0.0, 1.0) * z(17) * y(18) * z(19) + + double_complex(-1.0, 0.0) * z(17) * x(18) * z(19); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase7) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(11, 5, 4.0, 20); + cudaq::spin_op gold = + double_complex(0.0, 1.0) * z(3) * z(4) * x(5) * y(7) * z(9) * z(10) * + x(11) + + double_complex(-1.0, 0.0) * z(3) * y(5) * y(7) * z(9) * z(10) * x(11) + + double_complex(1.0, 0.0) * z(3) * z(4) * x(5) * y(7) * y(11) + + double_complex(0.0, 1.0) * z(3) * y(5) * y(7) * y(11); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase8) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(7, 9, 4.0, 20); + cudaq::spin_op gold = + double_complex(0.0, -1.0) * z(3) * z(5) * z(6) * y(7) * z(8) * x(9) * + x(11) + + double_complex(1.0, 0.0) * z(3) * z(5) * z(6) * y(7) * y(9) * x(11) + + double_complex(1.0, 0.0) * x(7) * z(8) * x(9) * x(11) + + double_complex(0.0, 1.0) * x(7) * y(9) * x(11); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase9) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(9, 15, 4.0, 20); + cudaq::spin_op gold = + double_complex(-1.0, 0.0) * y(9) * y(11) * z(13) * z(14) + + double_complex(0.0, -1.0) * z(8) * x(9) * y(11) * z(13) * z(14) + + double_complex(-1.0, 0.0) * z(7) * z(8) * x(9) * x(11) * z(15) + + double_complex(0.0, 1.0) * z(7) * y(9) * x(11) * z(15); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +} + +TEST(BravyiKitaev, testSRLCase10) { + using double_complex = std::complex; + using namespace cudaq::spin; + + auto result = cudaq::solvers::seeley_richard_love(3, 7, 4.0, 20); + cudaq::spin_op gold = + double_complex(0.0, -1.0) * z(1) * z(2) * y(3) * z(5) * z(6) + + double_complex(1.0, 0.0) * x(3) * z(5) * z(6) + + double_complex(-1.0, 0.0) * z(1) * z(2) * x(3) * z(7) + + double_complex(0.0, 1.0) * y(3) * z(7); + + auto [terms, residuals] = (result - gold).get_raw_data(); + for (auto r : residuals) + EXPECT_NEAR(std::abs(r), 0.0, 1e-4); +}