diff --git a/libs/solvers/lib/operators/molecule/fermion_compilers/jordan_wigner.cpp b/libs/solvers/lib/operators/molecule/fermion_compilers/jordan_wigner.cpp
index ed3dc16..c978e1d 100644
--- a/libs/solvers/lib/operators/molecule/fermion_compilers/jordan_wigner.cpp
+++ b/libs/solvers/lib/operators/molecule/fermion_compilers/jordan_wigner.cpp
@@ -5,6 +5,7 @@
  * This source code and the accompanying materials are made available under    *
  * the terms of the Apache License 2.0 which accompanies this distribution.    *
  ******************************************************************************/
+
 #include "cudaq/solvers/operators/molecule/fermion_compilers/jordan_wigner.h"
 
 #include <cppitertools/combinations.hpp>
@@ -16,426 +17,56 @@ using namespace cudaqx;
 
 namespace cudaq::solvers {
 
-cudaq::spin_op one_body(std::size_t p, std::size_t q,
-                        std::complex<double> coeff) {
-  if (p == q)
-    return 0.5 * coeff * cudaq::spin::i(p) - 0.5 * coeff * cudaq::spin::z(p);
-
-  if (p > q) {
-    std::swap(p, q);
-    coeff = std::conj(coeff);
-  }
-
-  std::vector<std::size_t> zIndices;
-  for (std::size_t i = p + 1; i < q; i++)
-    zIndices.push_back(i);
-
-  cudaq::spin_op parity = 1.;
-  for (auto i : zIndices)
-    parity *= cudaq::spin::z(i);
-
-  auto spin_hamiltonian =
-      0.5 * coeff.real() * cudaq::spin::x(p) * parity * cudaq::spin::x(q);
-  spin_hamiltonian +=
-      0.5 * coeff.real() * cudaq::spin::y(p) * parity * cudaq::spin::y(q);
-  spin_hamiltonian +=
-      0.5 * coeff.imag() * cudaq::spin::y(p) * parity * cudaq::spin::x(q);
-  spin_hamiltonian -=
-      0.5 * coeff.imag() * cudaq::spin::x(p) * parity * cudaq::spin::y(q);
-
-  return spin_hamiltonian;
-}
-
-cudaq::spin_op two_body(std::size_t p, std::size_t q, std::size_t r,
-                        std::size_t s, std::complex<double> coef) {
-  std::set<std::size_t> tmp{p, q, r, s};
-  if (tmp.size() == 2) {
-    auto spin_hamiltonian =
-        -0.25 * coef * cudaq::spin::i(p) * cudaq::spin::i(q);
-    if (p == r) {
-      spin_hamiltonian += 0.25 * coef * cudaq::spin::i(p) * cudaq::spin::z(q);
-      spin_hamiltonian += 0.25 * coef * cudaq::spin::z(p) * cudaq::spin::i(q);
-      spin_hamiltonian -= 0.25 * coef * cudaq::spin::z(p) * cudaq::spin::z(q);
-    } else if (q == r) {
-      spin_hamiltonian *= -1.;
-      spin_hamiltonian -= 0.25 * coef * cudaq::spin::i(p) * cudaq::spin::z(q);
-      spin_hamiltonian -= 0.25 * coef * cudaq::spin::z(p) * cudaq::spin::i(q);
-      spin_hamiltonian += 0.25 * coef * cudaq::spin::z(p) * cudaq::spin::z(q);
-    }
-    return spin_hamiltonian;
-  }
-
-  if (tmp.size() == 3) {
-    std::size_t a, b, c, d;
-    if (q == r) {
-      if (p > r) {
-        // a,b=s,p
-        a = s;
-        b = p;
-        coef = std::conj(coef);
-      } else {
-        // a,b=p,s
-        a = p;
-        b = s;
-      }
-      c = q;
-    } else if (q == s) {
-      if (p > r) {
-        // a,b=r,p
-        a = r;
-        b = p;
-        coef = -1.0 * std::conj(coef);
-      } else {
-        // a,b=p,r
-        a = p;
-        b = r;
-        coef *= -1.0;
-      }
-      c = q;
-    } else if (p == r) {
-      if (q > s) {
-        // a,b=s,q
-        a = s;
-        b = q;
-        coef = -1.0 * std::conj(coef);
-      } else {
-        // a,b=q,s
-        a = q;
-        b = s;
-        coef = -1.0 * coef;
-      }
-      c = p;
-    } else if (p == s) {
-      if (q > r) {
-        // a,b=r,q
-        a = r;
-        b = q;
-        coef = std::conj(coef);
-      } else {
-        // a,b=q,r
-        a = q;
-        b = r;
-      }
-      c = p;
-    }
-
-    std::vector<std::size_t> zIndices;
-    for (std::size_t i = a + 1; i < b; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parity = 1.;
-    for (auto i : zIndices)
-      parity *= cudaq::spin::z(i);
-
-    auto spin_hamiltonian = 0.25 * coef.real() * cudaq::spin::x(a) * parity *
-                            cudaq::spin::x(b) * cudaq::spin::i(c);
-    spin_hamiltonian += 0.25 * coef.real() * cudaq::spin::y(a) * parity *
-                        cudaq::spin::y(b) * cudaq::spin::i(c);
-    spin_hamiltonian += 0.25 * coef.imag() * cudaq::spin::y(a) * parity *
-                        cudaq::spin::x(b) * cudaq::spin::i(c);
-    spin_hamiltonian -= 0.25 * coef.imag() * cudaq::spin::x(a) * parity *
-                        cudaq::spin::y(b) * cudaq::spin::i(c);
-
-    spin_hamiltonian -= 0.25 * coef.real() * cudaq::spin::x(a) * parity *
-                        cudaq::spin::x(b) * cudaq::spin::z(c);
-    spin_hamiltonian -= 0.25 * coef.real() * cudaq::spin::y(a) * parity *
-                        cudaq::spin::y(b) * cudaq::spin::z(c);
-    spin_hamiltonian -= 0.25 * coef.imag() * cudaq::spin::y(a) * parity *
-                        cudaq::spin::x(b) * cudaq::spin::z(c);
-    spin_hamiltonian += 0.25 * coef.imag() * cudaq::spin::x(a) * parity *
-                        cudaq::spin::y(b) * cudaq::spin::z(c);
-    return spin_hamiltonian;
-  }
-
-  if ((p > q) ^ (r > s))
-    coef *= -1.0;
-
-  if (p < q && q < r && r < s) {
-    // a,b,c,d=p,q,r,s
-    auto a = p;
-    auto b = q;
-    auto c = r;
-    auto d = s;
-
-    std::vector<std::size_t> zIndices;
-    for (std::size_t i = a + 1; i < b; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parityA = 1.;
-    for (auto i : zIndices)
-      parityA *= cudaq::spin::z(i);
-
-    zIndices.clear();
-    for (std::size_t i = c + 1; i < d; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parityB = 1.;
-    for (auto i : zIndices)
-      parityB *= cudaq::spin::z(i);
-
-    auto spin_hamiltonian = -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                            cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                            cudaq::spin::x(d);
-    spin_hamiltonian -= -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian -= -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    return spin_hamiltonian;
-  }
-
-  if (p < r && r < q && q < s) {
-    // a,b,c,d=p,r,q,s
-    auto a = p;
-    auto b = r;
-    auto c = q;
-    auto d = s;
-
-    std::vector<std::size_t> zIndices;
-    for (std::size_t i = a + 1; i < b; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parityA = 1.;
-    for (auto i : zIndices)
-      parityA *= cudaq::spin::z(i);
-
-    zIndices.clear();
-    for (std::size_t i = c + 1; i < d; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parityB = 1.;
-    for (auto i : zIndices)
-      parityB *= cudaq::spin::z(i);
-
-    auto spin_hamiltonian = -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                            cudaq::spin::x(b) * cudaq::spin::x(c) * parityA *
-                            cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::y(d);
-
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::x(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityA *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityA *
-                        cudaq::spin::x(d);
-    return spin_hamiltonian;
-  }
-
-  if (p < r && r < s && s < q) {
-    // a,b,c,d=p,r,s,q
-    auto a = p;
-    auto b = r;
-    auto c = s;
-    auto d = q;
-
-    std::vector<std::size_t> zIndices;
-    for (std::size_t i = a + 1; i < b; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parityA = 1.;
-    for (auto i : zIndices)
-      parityA *= cudaq::spin::z(i);
-
-    zIndices.clear();
-    for (std::size_t i = c + 1; i < d; i++)
-      zIndices.push_back(i);
-
-    cudaq::spin_op parityB = 1.;
-    for (auto i : zIndices)
-      parityB *= cudaq::spin::z(i);
-
-    auto spin_hamiltonian = -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                            cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                            cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= -0.125 * coef.real() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian -= -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += -0.125 * coef.real() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::x(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::x(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::x(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian -= 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::x(c) * parityB *
-                        cudaq::spin::y(d);
-    spin_hamiltonian += 0.125 * coef.imag() * cudaq::spin::y(a) * parityA *
-                        cudaq::spin::y(b) * cudaq::spin::y(c) * parityB *
-                        cudaq::spin::x(d);
-
-    return spin_hamiltonian;
-  }
-
-  throw std::runtime_error(
-      "Invalid condition in two_body jordan wigner function.");
-}
-
 cudaq::spin_op jordan_wigner::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 spin_hamiltonian = constant * cudaq::spin_op();
   std::size_t nqubit = hpq.shape()[0];
-
-  double tolerance =
+  double tol =
       options.get<double>(std::vector<std::string>{"tolerance", "tol"}, 1e-15);
 
-  for (auto p : cudaq::range(nqubit)) {
-    auto coef = hpq.at({p, p});
-    if (std::fabs(coef) > tolerance)
-      spin_hamiltonian += one_body(p, p, coef);
-  }
-
-  std::vector<std::vector<std::size_t>> next;
-  for (auto &&combo : iter::combinations(cudaq::range(nqubit), 2)) {
-    auto p = combo[0];
-    auto q = combo[1];
-    next.push_back({p, q});
-    auto coef = 0.5 * (hpq.at({p, q}) + std::conj(hpq.at({q, p})));
-    if (std::fabs(coef) > tolerance)
-      spin_hamiltonian += one_body(p, q, coef);
-
-    coef = hpqrs.at({p, q, p, q}) + hpqrs.at({q, p, q, p});
-    if (std::fabs(coef) > tolerance)
-      spin_hamiltonian += two_body(p, q, p, q, coef);
-
-    coef = hpqrs.at({p, q, q, p}) + hpqrs.at({q, p, p, q});
-    if (std::fabs(coef) > tolerance)
-      spin_hamiltonian += two_body(p, q, q, p, coef);
-  }
-
-  for (auto combo : iter::combinations(next, 2)) {
-    auto p = combo[0][0];
-    auto q = combo[0][1];
-    auto r = combo[1][0];
-    auto s = combo[1][1];
-    auto coef =
-        0.5 * (hpqrs.at({p, q, r, s}) + std::conj(hpqrs.at({s, r, q, p})) -
-               hpqrs.at({q, p, r, s}) - std::conj(hpqrs.at({s, r, p, q})) -
-               hpqrs.at({p, q, s, r}) - std::conj(hpqrs.at({r, s, q, p})) +
-               hpqrs.at({q, p, s, r}) + std::conj(hpqrs.at({r, s, p, q})));
-
-    if (std::fabs(coef) > tolerance)
-      spin_hamiltonian += two_body(p, q, r, s, coef);
-  }
+  auto is_complex_zero = [tol](const std::complex<double> &z) {
+    return std::abs(z.real()) < tol && std::abs(z.imag()) < tol;
+  };
+
+  auto adag = [](std::size_t numQubits, std::size_t j) {
+    cudaq::spin_op zprod(numQubits);
+    for (std::size_t k = 0; k < j; k++)
+      zprod *= cudaq::spin::z(k);
+    return 0.5 * zprod *
+           (cudaq::spin::x(j) - std::complex<double>{0, 1} * cudaq::spin::y(j));
+  };
+
+  auto a = [](std::size_t numQubits, std::size_t j) {
+    cudaq::spin_op zprod(numQubits);
+    for (std::size_t k = 0; k < j; k++)
+      zprod *= cudaq::spin::z(k);
+    return 0.5 * zprod *
+           (cudaq::spin::x(j) + std::complex<double>{0, 1} * cudaq::spin::y(j));
+  };
+
+  for (std::size_t i = 0; i < hpq.shape()[0]; i++)
+    for (std::size_t j = 0; j < hpq.shape()[1]; j++)
+      if (!is_complex_zero(hpq.at({i, j})))
+        spin_hamiltonian += hpq.at({i, j}) * adag(nqubit, i) * a(nqubit, j);
+
+  for (std::size_t i = 0; i < hpqrs.shape()[0]; i++)
+    for (std::size_t j = 0; j < hpqrs.shape()[1]; j++)
+      for (std::size_t k = 0; k < hpqrs.shape()[0]; k++)
+        for (std::size_t l = 0; l < hpqrs.shape()[1]; l++)
+          if (!is_complex_zero(hpqrs.at({i, j, k, l})))
+            spin_hamiltonian += hpqrs.at({i, j, k, l}) * adag(nqubit, i) *
+                                adag(nqubit, j) * a(nqubit, k) * a(nqubit, l);
 
   // Remove terms with 0.0 coefficient
   std::vector<cudaq::spin_op> nonZeros;
-  for (auto term : spin_hamiltonian) {
+  for (auto &term : spin_hamiltonian) {
     auto coeff = term.get_coefficient();
-    if (std::fabs(coeff) > tolerance)
+    if (std::fabs(coeff) > tol)
       nonZeros.push_back(term);
   }
   auto op = nonZeros[0];
@@ -444,4 +75,4 @@ cudaq::spin_op jordan_wigner::generate(const double constant,
 
   return op;
 }
-} // namespace cudaq::solvers
\ No newline at end of file
+} // namespace cudaq::solvers
diff --git a/libs/solvers/python/tests/test_jordan_wigner.py b/libs/solvers/python/tests/test_jordan_wigner.py
new file mode 100644
index 0000000..4037907
--- /dev/null
+++ b/libs/solvers/python/tests/test_jordan_wigner.py
@@ -0,0 +1,115 @@
+# ============================================================================ #
+# Copyright (c) 2025 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.     #
+# ============================================================================ #
+
+import os
+import pytest
+import cudaq, cudaq_solvers as solvers
+from pyscf import gto, scf, fci
+import numpy as np
+
+
+def extract_words(hamiltonian: cudaq.SpinOperator):
+    result = []
+    hamiltonian.for_each_term(lambda term: result.append(term.to_string(False)))
+    return result
+
+
+def extract_coefficients(hamiltonian: cudaq.SpinOperator):
+    result = []
+    hamiltonian.for_each_term(
+        lambda term: result.append(term.get_coefficient()))
+    return result
+
+
+def extract_spin_op_to_dict(op: cudaq.SpinOperator) -> dict:
+    d = {}
+    coeffs = extract_coefficients(op)
+    words = extract_words(op)
+    for c, w in zip(coeffs, words):
+        d[w] = c
+    return d
+
+
+def jw_molecule_compare_hamiltonians_test(xyz):
+    # Compute energy using CUDA-Q/OpenFermion
+    of_hamiltonian, data = cudaq.chemistry.create_molecular_hamiltonian(
+        xyz, 'sto-3g', 1, 0)
+
+    # Compute energy using CUDA-QX
+    molecule = solvers.create_molecule(xyz, 'sto-3g', 0, 0, casci=True)
+    cqx_op = solvers.jordan_wigner(
+        molecule.hpq,
+        molecule.hpqrs,
+        core_energy=molecule.energies['nuclear_energy'],
+        tol=1e-12)
+
+    of_hamiltonian_dict = extract_spin_op_to_dict(of_hamiltonian)
+    cqx_op_dict = extract_spin_op_to_dict(cqx_op)
+
+    for k in of_hamiltonian_dict.keys():
+        assert (k in cqx_op_dict.keys())
+    for k in cqx_op_dict.keys():
+        assert (k in of_hamiltonian_dict.keys())
+
+    for k in of_hamiltonian_dict.keys():
+        np.isclose(of_hamiltonian_dict[k], cqx_op_dict[k], 1e-12)
+
+
+def jw_molecule_test(xyz):
+    # Compute FCI energy
+    mol = gto.M(atom=xyz, basis='sto-3g', symmetry=False)
+    mf = scf.RHF(mol).run()
+    fci_energy = fci.FCI(mf).kernel()[0]
+    print(f'FCI energy:            {fci_energy}')
+
+    # Compute energy using CUDA-Q/OpenFermion
+    of_hamiltonian, data = cudaq.chemistry.create_molecular_hamiltonian(
+        xyz, 'sto-3g', 1, 0)
+    of_energy = np.min(np.linalg.eigvals(of_hamiltonian.to_matrix()))
+    print(f'OpenFermion energy:    {of_energy.real}')
+
+    # Compute energy using CUDA-QX
+    molecule = solvers.create_molecule(xyz, 'sto-3g', 0, 0, casci=True)
+    op = solvers.jordan_wigner(molecule.hpq,
+                               molecule.hpqrs,
+                               core_energy=molecule.energies['nuclear_energy'])
+    assert op == molecule.hamiltonian
+    assert of_hamiltonian == molecule.hamiltonian
+
+    # FIXME - why do we need to call eigvals again if we can just assert the equality checks above?
+    cudaqx_eig = np.min(np.linalg.eigvals(op.to_matrix()))
+    print(f'CUDA-QX energy:        {cudaqx_eig.real}')
+    assert np.isclose(cudaqx_eig, of_energy.real, atol=1e-4)
+
+    num_terms = of_hamiltonian.get_term_count()
+    print(f'Number of terms in CUDA-Q/OpenFermion: {num_terms}')
+    num_terms = molecule.hamiltonian.get_term_count()
+    print(f'Number of terms in CUDA-QX           : {num_terms}')
+
+
+def test_ground_state():
+    xyz = [('H', (0., 0., 0.)), ('H', (0., 0., .7474))]
+    jw_molecule_compare_hamiltonians_test(xyz)
+    jw_molecule_test(xyz)
+    xyz = [('H', (0., 0., 0.)), ('H', (0., 0., .7474)), ('H', (1., 0., 0.)),
+           ('H', (1., 0., .7474))]
+    jw_molecule_compare_hamiltonians_test(xyz)
+    jw_molecule_test(xyz)
+    xyz = [('H', (0., 0., 0.)), ('H', (1.0, 0., 0.)),
+           ('H', (0.322, 2.592, 0.1)), ('H', (1.2825, 2.292, 0.1))]
+    jw_molecule_compare_hamiltonians_test(xyz)
+    jw_molecule_test(xyz)
+    xyz = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.1774))]
+    jw_molecule_compare_hamiltonians_test(xyz)
+    jw_molecule_test(xyz)
+    xyz = [('O', (0.000000, 0.000000, 0.000000)),
+           ('H', (0.757000, 0.586000, 0.000000)),
+           ('H', (-0.757000, 0.586000, 0.000000))]
+    jw_molecule_compare_hamiltonians_test(xyz)
+    # This is commented out by default because it is a very long test.
+    # jw_molecule_test(xyz)