wip

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

cudaqlib/operators.h

@@ -9,4 +9,5 @@
 
 #include "operators/chemistry/MoleculePackageDriver.h"
 #include "operators/fermion/fermion_compiler.h"
+#include "operators/fermion/fermion_operators.h"
 #include "operators/fermion/transformations/jordan_wigner.h"

cudaqlib/operators/chemistry/molecule.cpp

@@ -111,7 +111,31 @@ void molecule_options::dump() {
 
 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);
+  if (p == q)
+    return 0.5 * spin::i(numQubits - 1) * 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 spin::i(numQubits - 1) * ret;
 }
 
 } // namespace cudaq::operators

cudaqlib/operators/fermion/transformations/jordan_wigner.cpp

@@ -383,13 +383,14 @@ spin_op jordan_wigner::generate(const fermion_op &fermionOp) {
 }
 
 spin_op jordan_wigner::adag(std::size_t numQubits, std::size_t p) {
-  spin_op zprod(numQubits);
+  spin_op zprod; // = spin::i(numQubits - 1);
   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);
+  spin_op zprod; // = spin::i(numQubits - 1);
   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));

python/bindings/operators/py_operators.cpp

@@ -80,7 +80,12 @@ void bindOperators(py::module &mod) {
       .def_readonly("n_electrons", &molecular_hamiltonian::n_electrons)
       .def_readonly("n_orbitals", &molecular_hamiltonian::n_orbitals);
 
-  operators.def("one_particle_op", &one_particle_op, "");
+  operators.def(
+      "one_particle_op",
+      [](std::size_t numQ, std::size_t p, std::size_t q) {
+        return operators::one_particle_op(numQ, p, q);
+      },
+      "");
 
   operators.def(
       "create_molecule",

tests/operators/OperatorsTester.cpp

@@ -505,11 +505,7 @@ H -0.4691 -0.7570 0.0
 }
 
 TEST(OperatorsTester, checkOneParticleOp) {
-  std::size_t numQubits = 4;
-
   auto op = cudaq::operators::one_particle_op(4, 1, 1);
-  op.dump();
-
-  op = cudaq::operators::one_particle_op(5, 2, 4);
+  EXPECT_EQ(op.num_terms(), 2);
   op.dump();
 }