Better one_particle_op implemenation

Signed-off-by: Alex McCaskey <[email protected]>

cudaqlib/gse/adapt/adapt.h

@@ -74,7 +74,7 @@ auto adapt_vqe(const InitialState &initialState, const spin_op &H,
     auto maxOpIdx = std::distance(gradients.begin(), iter);
 
     // Convergence is reached if gradient values are small
-    if (std::sqrt(std::fabs(norm)) < options.grad_norm_magnitude ||
+    if (std::sqrt(std::fabs(norm)) < options.grad_norm_tolerance ||
         std::fabs(lastNorm - norm) < options.grad_norm_tolerance)
       break;
 

cudaqlib/gse/utils/operator_pool.h

@@ -52,20 +52,4 @@ CUDAQ_DEFINE_EXTENSION_IMPL(operator_pool)
     return std::make_unique<TYPE>();                                           \
   }
 
-namespace fermion {
-inline spin_op adag(std::size_t numQubits, std::size_t j) {
-  spin_op zprod(numQubits);
-  for (std::size_t k = 0; k < j; k++)
-    zprod *= spin::z(k);
-  return 0.5 * zprod * (spin::x(j) - std::complex<double>{0, 1} * spin::y(j));
-}
-
-inline spin_op a(std::size_t numQubits, std::size_t j) {
-  spin_op zprod(numQubits);
-  for (std::size_t k = 0; k < j; k++)
-    zprod *= spin::z(k);
-  return 0.5 * zprod * (spin::x(j) + std::complex<double>{0, 1} * spin::y(j));
-}
-} // namespace fermion
-
 } // namespace cudaq

cudaqlib/gse/utils/pools/CMakeLists.txt

@@ -12,7 +12,7 @@ add_library(${LIBRARY_NAME} SHARED uccsd_pool.cpp spin_complement_gsd.cpp)
 target_include_directories(
   ${LIBRARY_NAME}
   PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}>)
-target_link_libraries(${LIBRARY_NAME} PRIVATE cudaq::cudaq-spin)
+target_link_libraries(${LIBRARY_NAME} PRIVATE cudaq::cudaq-spin cudaq-operators)
 
 install(
     TARGETS ${LIBRARY_NAME}

cudaqlib/gse/utils/pools/spin_complement_gsd.cpp

@@ -7,6 +7,7 @@
  ******************************************************************************/
 
 #include "spin_complement_gsd.h"
+#include "cudaqlib/operators/fermion/fermion_operators.h"
 
 namespace cudaq {
 
@@ -32,11 +33,11 @@ std::vector<spin_op> spin_complement_gsd::generate(
       if (p >= q)
         continue;
       auto oneElectron =
-          fermion::adag(numQubits, q) * fermion::a(numQubits, p) -
-          fermion::adag(numQubits, p) * fermion::a(numQubits, q);
+          operators::adag(numQubits, q) * operators::a(numQubits, p) -
+          operators::adag(numQubits, p) * operators::a(numQubits, q);
       oneElectron +=
-          fermion::adag(numQubits, q + 1) * fermion::a(numQubits, p + 1) -
-          fermion::adag(numQubits, p + 1) * fermion::a(numQubits, q + 1);
+          operators::adag(numQubits, q + 1) * operators::a(numQubits, p + 1) -
+          operators::adag(numQubits, p + 1) * operators::a(numQubits, q + 1);
 
       std::unordered_map<spin_op::spin_op_term, std::complex<double>> terms;
       oneElectron.for_each_term([&](spin_op &term) {
@@ -69,17 +70,18 @@ std::vector<spin_op> spin_complement_gsd::generate(
             continue;
 
           auto twoElectron =
-              fermion::adag(numQubits, r) * fermion::a(numQubits, p) *
-                  fermion::adag(numQubits, s) * fermion::a(numQubits, q) -
-              fermion::adag(numQubits, q) * fermion::a(numQubits, s) *
-                  fermion::adag(numQubits, p) * fermion::a(numQubits, r);
-          twoElectron +=
-              fermion::adag(numQubits, r + 1) * fermion::a(numQubits, p + 1) *
-                  fermion::adag(numQubits, s + 1) *
-                  fermion::a(numQubits, q + 1) -
-              fermion::adag(numQubits, q + 1) * fermion::a(numQubits, s + 1) *
-                  fermion::adag(numQubits, p + 1) *
-                  fermion::a(numQubits, r + 1);
+              operators::adag(numQubits, r) * operators::a(numQubits, p) *
+                  operators::adag(numQubits, s) * operators::a(numQubits, q) -
+              operators::adag(numQubits, q) * operators::a(numQubits, s) *
+                  operators::adag(numQubits, p) * operators::a(numQubits, r);
+          twoElectron += operators::adag(numQubits, r + 1) *
+                             operators::a(numQubits, p + 1) *
+                             operators::adag(numQubits, s + 1) *
+                             operators::a(numQubits, q + 1) -
+                         operators::adag(numQubits, q + 1) *
+                             operators::a(numQubits, s + 1) *
+                             operators::adag(numQubits, p + 1) *
+                             operators::a(numQubits, r + 1);
 
           std::unordered_map<spin_op::spin_op_term, std::complex<double>> terms;
           twoElectron.for_each_term([&](spin_op &term) {
@@ -114,20 +116,21 @@ std::vector<spin_op> spin_complement_gsd::generate(
             continue;
 
           auto twoElectron =
-              fermion::adag(numQubits, r) * fermion::a(numQubits, p) *
-                  fermion::adag(numQubits, s) * fermion::a(numQubits, q) -
-              fermion::adag(numQubits, q) * fermion::a(numQubits, s) *
-                  fermion::adag(numQubits, p) * fermion::a(numQubits, r);
+              operators::adag(numQubits, r) * operators::a(numQubits, p) *
+                  operators::adag(numQubits, s) * operators::a(numQubits, q) -
+              operators::adag(numQubits, q) * operators::a(numQubits, s) *
+                  operators::adag(numQubits, p) * operators::a(numQubits, r);
           if (p > q)
             continue;
 
-          twoElectron +=
-              fermion::adag(numQubits, s - 1) * fermion::a(numQubits, q - 1) *
-                  fermion::adag(numQubits, r + 1) *
-                  fermion::a(numQubits, p + 1) -
-              fermion::adag(numQubits, p + 1) * fermion::a(numQubits, r + 1) *
-                  fermion::adag(numQubits, q - 1) *
-                  fermion::a(numQubits, s - 1);
+          twoElectron += operators::adag(numQubits, s - 1) *
+                             operators::a(numQubits, q - 1) *
+                             operators::adag(numQubits, r + 1) *
+                             operators::a(numQubits, p + 1) -
+                         operators::adag(numQubits, p + 1) *
+                             operators::a(numQubits, r + 1) *
+                             operators::adag(numQubits, q - 1) *
+                             operators::a(numQubits, s - 1);
           std::unordered_map<spin_op::spin_op_term, std::complex<double>> terms;
           twoElectron.for_each_term([&](spin_op &term) {
             auto coeff = term.get_coefficient();

cudaqlib/operators.h

@@ -8,5 +8,5 @@
 #pragma once
 
 #include "operators/chemistry/MoleculePackageDriver.h"
-#include "operators/fermion/fermion_to_spin.h"
+#include "operators/fermion/fermion_compiler.h"
 #include "operators/fermion/transformations/jordan_wigner.h"

cudaqlib/operators/CMakeLists.txt

@@ -12,6 +12,7 @@ add_library(${LIBRARY_NAME} SHARED
   # Base implementation files
   chemistry/molecule.cpp
   fermion/fermion_op.cpp
+  fermion/fermion_operators.cpp
   # fermion transformations
   fermion/transformations/jordan_wigner.cpp
   # chemistry package drivers

cudaqlib/operators/chemistry/drivers/pyscf/ExternalPySCFDriver.cpp

@@ -9,7 +9,7 @@
 #include "nlohmann/json.hpp"
 
 #include "cudaqlib/operators/chemistry/MoleculePackageDriver.h"
-#include "cudaqlib/operators/fermion/fermion_to_spin.h"
+#include "cudaqlib/operators/fermion/fermion_compiler.h"
 #include "cudaqlib/utils/library_utils.h"
 
 #include <filesystem>
@@ -118,7 +118,7 @@ class external_pyscf
     fermionOp.hpqrs.set_data(hpqrsValues);
 
     // Transform to a spin operator
-    auto transform = fermion_to_spin::get(options.fermion_to_string);
+    auto transform = fermion_compiler::get(options.fermion_to_string);
     auto spinHamiltonian = transform->generate(fermionOp);
 
     // Return the molecular hamiltonian

cudaqlib/operators/chemistry/drivers/pyscf/RESTPySCFDriver.cpp

@@ -9,7 +9,7 @@
 #include "nlohmann/json.hpp"
 
 #include "cudaqlib/operators/chemistry/MoleculePackageDriver.h"
-#include "cudaqlib/operators/fermion/fermion_to_spin.h"
+#include "cudaqlib/operators/fermion/fermion_compiler.h"
 #include "cudaqlib/utils/library_utils.h"
 
 #include "common/RestClient.h"
@@ -109,7 +109,7 @@ class RESTPySCFDriver
     fermionOp.hpqrs.set_data(hpqrsValues);
 
     // Transform to a spin operator
-    auto transform = fermion_to_spin::get(options.fermion_to_string);
+    auto transform = fermion_compiler::get(options.fermion_to_string);
     auto spinHamiltonian = transform->generate(fermionOp);
 
     // Return the molecular hamiltonian

cudaqlib/operators/chemistry/molecule.cpp

@@ -7,6 +7,7 @@
  ******************************************************************************/
 #include "molecule.h"
 #include "MoleculePackageDriver.h"
+#include "cudaqlib/operators/fermion/fermion_operators.h"
 
 #include <fstream>
 #include <iostream>
@@ -108,33 +109,9 @@ void molecule_options::dump() {
             << "\n]\n";
 }
 
-spin_op one_particle_op(std::size_t p, std::size_t q) {
-
-  if (p == q)
-    return 0.5 * spin::i(p) - 0.5 * spin::z(p);
-
-  std::complex<double> coeff(0., 1.);
-  double m = -.25;
-  if (p > q) {
-    std::swap(p, q);
-    coeff = std::conj(coeff);
-  }
-
-  std::vector<std::size_t> z_indices;
-  for (auto i : cudaq::range((long)p + 1, (long)q))
-    z_indices.push_back(i);
-
-  auto parity = spin::z(z_indices.front());
-  for (std::size_t i = 1; i < z_indices.size(); i++) {
-    parity *= spin::z(i);
-  }
-
-  auto ret = m * spin::x(p) * parity * spin::x(q);
-
-  ret += m * spin::y(p) * parity * spin::y(q);
-  ret -= coeff * m * spin::y(p) * parity * spin::x(q);
-  ret += coeff * m * spin::x(p) * parity * spin::y(q);
-  return ret;
+spin_op one_particle_op(std::size_t numQubits, std::size_t p, std::size_t q,
+                        const std::string fermionCompiler) {
+  return adag(numQubits, p, fermionCompiler) * a(numQubits, q, fermionCompiler);
 }
 
-} // namespace cudaq::operators
+} // namespace cudaq::operators

cudaqlib/operators/chemistry/molecule.h

@@ -83,6 +83,7 @@ create_molecule(const molecular_geometry &geometry, const std::string &basis,
                 int spin, int charge,
                 molecule_options options = molecule_options());
 
-spin_op one_particle_op(std::size_t p, std::size_t q);
+spin_op one_particle_op(std::size_t numQubits, std::size_t p, std::size_t q,
+                        const std::string fermionCompiler = "jordan_wigner");
 
 } // namespace cudaq::operators

cudaqlib/operators/fermion/fermion_to_spin.h → cudaqlib/operators/fermion/fermion_compiler.h

@@ -1,5 +1,5 @@
 /*******************************************************************************
- * Copyright (c) 2022 - 2023 NVIDIA Corporation & Affiliates.                  *
+ * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates.                  *
  * All rights reserved.                                                        *
  *                                                                             *
  * This source code and the accompanying materials are made available under    *
@@ -14,22 +14,24 @@
 
 namespace cudaq::operators {
 
-/// @brief The `fermion_to_spin` type serves as a base class defining
+/// @brief The `fermion_compiler` type serves as a base class defining
 /// and interface for clients to map `fermion_op` instances to
 /// `cudaq::spin_op` instances.
-class fermion_to_spin : public extension_point<fermion_to_spin> {
+class fermion_compiler : public extension_point<fermion_compiler> {
 public:
   /// @brief Given a fermionic representation of an operator
   /// generate an equivalent operator on spins.
   virtual spin_op generate(const fermion_op &fermionOp) = 0;
-  virtual ~fermion_to_spin() {}
+  virtual spin_op adag(std::size_t numQubits, std::size_t p) = 0;
+  virtual spin_op a(std::size_t numQubits, std::size_t q) = 0;
+  virtual ~fermion_compiler() {}
 };
 
-CUDAQ_DEFINE_EXTENSION_IMPL(fermion_to_spin)
+CUDAQ_DEFINE_EXTENSION_IMPL(fermion_compiler)
 
 #define CUDAQ_REGISTER_FERMION_TO_SPIN(TYPE)                                   \
   static inline const std::string class_identifier = #TYPE;                    \
-  static std::unique_ptr<fermion_to_spin> create() {                           \
+  static std::unique_ptr<fermion_compiler> create() {                          \
     return std::make_unique<TYPE>();                                           \
   }
 

cudaqlib/operators/fermion/fermion_operators.cpp

@@ -0,0 +1,22 @@
+/*******************************************************************************
+ * Copyright (c) 2022 - 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 "fermion_operators.h"
+#include "fermion_compiler.h"
+
+namespace cudaq::operators {
+spin_op adag(std::size_t numQubits, std::size_t p,
+             const std::string transform) {
+  auto transformer = fermion_compiler::get(transform);
+  return transformer->adag(numQubits, p);
+}
+spin_op a(std::size_t numQubits, std::size_t q, const std::string transform) {
+  auto transformer = fermion_compiler::get(transform);
+  return transformer->adag(numQubits, q);
+}
+
+} // namespace cudaq::operators

cudaqlib/operators/fermion/fermion_operators.h

@@ -0,0 +1,17 @@
+/*******************************************************************************
+ * Copyright (c) 2022 - 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/spin_op.h"
+
+namespace cudaq::operators {
+spin_op adag(std::size_t numQubits, std::size_t p,
+             const std::string transform = "jordan_wigner");
+spin_op a(std::size_t numQubits, std::size_t q,
+          const std::string transform = "jordan_wigner");
+} // namespace cudaq::operators

cudaqlib/operators/fermion/transformations/jordan_wigner.cpp

@@ -381,4 +381,17 @@ spin_op jordan_wigner::generate(const fermion_op &fermionOp) {
 
   return op;
 }
+
+spin_op jordan_wigner::adag(std::size_t numQubits, std::size_t p) {
+  spin_op zprod(numQubits);
+  for (std::size_t k = 0; k < p; k++)
+    zprod *= spin::z(k);
+  return 0.5 * zprod * (spin::x(p) - std::complex<double>{0, 1} * spin::y(p));
+}
+spin_op jordan_wigner::a(std::size_t numQubits, std::size_t q) {
+  spin_op zprod(numQubits);
+  for (std::size_t k = 0; k < q; k++)
+    zprod *= spin::z(k);
+  return 0.5 * zprod * (spin::x(q) + std::complex<double>{0, 1} * spin::y(q));
+}
 } // namespace cudaq::operators

cudaqlib/operators/fermion/transformations/jordan_wigner.h

@@ -7,15 +7,17 @@
  ******************************************************************************/
 #pragma once
 
-#include "../fermion_to_spin.h"
+#include "../fermion_compiler.h"
 
 namespace cudaq::operators {
 
 /// @brief Map fermionic operators to spin operators via the
 /// Jordan-Wigner transformation.
-class jordan_wigner : public details::fermion_to_spin_impl<jordan_wigner> {
+class jordan_wigner : public details::fermion_compiler_impl<jordan_wigner> {
 public:
   spin_op generate(const fermion_op &fermionOp) override;
+  spin_op adag(std::size_t numQubits, std::size_t p) override;
+  spin_op a(std::size_t numQubits, std::size_t q) override;
   CUDAQ_REGISTER_FERMION_TO_SPIN(jordan_wigner)
 };
 } // namespace cudaq::operators

python/bindings/operators/py_operators.cpp

@@ -23,14 +23,14 @@ void bindOperators(py::module &mod) {
   auto operators = mod.def_submodule("operators");
 
   operators.def("jordan_wigner", [](fermion_op &op) {
-    return fermion_to_spin::get("jordan_wigner")->generate(op);
+    return fermion_compiler::get("jordan_wigner")->generate(op);
   });
 
-  py::class_<fermion_to_spin>(operators, "fermion_to_spin")
+  py::class_<fermion_compiler>(operators, "fermion_compiler")
       .def_static(
           "get",
-          [](const std::string &name) { return fermion_to_spin::get(name); })
-      .def("generate", &fermion_to_spin::generate, "");
+          [](const std::string &name) { return fermion_compiler::get(name); })
+      .def("generate", &fermion_compiler::generate, "");
 
   py::class_<one_body_integrals>(operators, "OneBodyIntegrals",
                                  py::buffer_protocol())

python/tests/CMakeLists.txt

@@ -0,0 +1,8 @@
+# ============================================================================ #
+# Copyright (c) 2022 - 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.     #
+# ============================================================================ #
+message(STATUS "STATUS WE IN HERE")

tests/gse/AdaptTester.cpp

@@ -27,7 +27,6 @@ TEST(GSETester, checkSimpleAdapt) {
 
   auto pool = cudaq::operator_pool::get("spin_complement_gsd");
   auto poolList = pool->generate({{"num-orbitals", h.num_qubits() / 2}});
-
   auto initialState = [&](cudaq::qvector<> &q) __qpu__ {
     for (std::size_t i = 0; i < 2; i++)
       x(q[i]);