Implements fast handling of Pauli string

CMakeLists.txt

@@ -44,14 +44,17 @@ include_directories(${SOURCE_DIR})
 add_subdirectory(lib/pybind11)
 add_subdirectory(lib/fmt)
 pybind11_add_module(qforte ${SOURCES} "${SOURCE_DIR}/bindings.cc"
+    "${SOURCE_DIR}/bitstring.cc"
     "${SOURCE_DIR}/helpers.cc"
     "${SOURCE_DIR}/make_gate.cc"
-    "${SOURCE_DIR}/qubit_basis.cc"
     "${SOURCE_DIR}/circuit.cc"
     "${SOURCE_DIR}/computer.cc"
     "${SOURCE_DIR}/gate.cc"
+    "${SOURCE_DIR}/pauli_string.cc"
     "${SOURCE_DIR}/qubit_operator.cc"
     "${SOURCE_DIR}/qubit_op_pool.cc"
+    "${SOURCE_DIR}/pauli_string_vector.cc"
+    "${SOURCE_DIR}/qubit_basis.cc"
     "${SOURCE_DIR}/sq_operator.cc"
     "${SOURCE_DIR}/sq_op_pool.cc"
     "${SOURCE_DIR}/sparse_tensor.cc"
@@ -83,15 +86,18 @@ include_directories(lib/fmt/include)
 
 add_executable("${PROJECT_NAME}_benchmarks"
     benchmarks/benchmarks.cc
+    "${SOURCE_DIR}/bitstring.cc"
     "${SOURCE_DIR}/helpers.cc"
     "${SOURCE_DIR}/make_gate.cc"
-    "${SOURCE_DIR}/qubit_basis.cc"
     "${SOURCE_DIR}/circuit.cc"
     "${SOURCE_DIR}/computer.cc"
     "${SOURCE_DIR}/gate.cc"
+    "${SOURCE_DIR}/pauli_string.cc"
     "${SOURCE_DIR}/qubit_operator.cc"
     "${SOURCE_DIR}/qubit_op_pool.cc"
+    "${SOURCE_DIR}/pauli_string_vector.cc"
     "${SOURCE_DIR}/sparse_tensor.cc"
+    "${SOURCE_DIR}/qubit_basis.cc"
     "${SOURCE_DIR}/timer.cc")
 
 target_link_libraries(qforte_benchmarks PRIVATE fmt-header-only)

src/qforte/abc/algorithm.py

@@ -30,6 +30,11 @@ class Algorithm(ABC):
         measurment (unphysical for quantum computer). Most algorithms only
         have a fast implentation.
 
+    _fast_Pauli : bool
+        If true, the code uses the bitstring representation of Pauli strings
+        and qubit basis states in the computationally demanding parts of the
+        algorithm.
+
     _trotter_order : int
         The Trotter order to use for exponentiated operators.
         (exact in the infinite limit).
@@ -71,6 +76,7 @@ def __init__(self,
                  trotter_order=1,
                  trotter_number=1,
                  fast=True,
+                 fast_Pauli=True,
                  verbose=False,
                  print_summary_file=False,
                  **kwargs):
@@ -116,6 +122,7 @@ def __init__(self,
         self._trotter_order = trotter_order
         self._trotter_number = trotter_number
         self._fast = fast
+        self._fast_Pauli = fast_Pauli
         self._verbose = verbose
         self._print_summary_file = print_summary_file
 
@@ -270,7 +277,7 @@ def fill_pool(self):
                     temp_sq_pool.add(sq_operator[0], sq_operator[1])
             self._pool_obj = temp_sq_pool
 
-        self._Nm = [len(operator.jw_transform().terms()) for _, operator in self._pool_obj]
+        self._Nm = [len(operator.jw_transform(self._fast_Pauli).terms()) for _, operator in self._pool_obj]
 
     def measure_energy(self, Ucirc):
         """

src/qforte/abc/ansatz.py

@@ -31,7 +31,7 @@ def ansatz_circuit(self, amplitudes=None):
         for tamp, top in zip(tamps, self._tops):
             temp_pool.add(tamp, self._pool_obj[top][1])
 
-        A = temp_pool.get_qubit_operator('commuting_grp_lex')
+        A = temp_pool.get_qubit_operator('commuting_grp_lex', True, self._fast_Pauli)
 
         U, phase1 = trotterize(A, trotter_number=self._trotter_number)
         if phase1 != 1.0 + 0.0j:

src/qforte/abc/uccvqeabc.py

@@ -78,7 +78,7 @@ def get_num_commut_measurements(self):
 
     def fill_commutator_pool(self):
         print('\n\n==> Building commutator pool for gradient measurement.')
-        self._commutator_pool = self._pool_obj.get_qubit_op_pool()
+        self._commutator_pool = self._pool_obj.get_qubit_op_pool(self._fast_Pauli)
         self._commutator_pool.join_as_commutator(self._qb_ham)
         print('==> Commutator pool construction complete.')
 
@@ -145,7 +145,7 @@ def measure_gradient(self, params=None):
         mu = M-1
 
         # find <sing_N | K_N | psi_N>
-        Kmu_prev = self._pool_obj[self._tops[mu]][1].jw_transform()
+        Kmu_prev = self._pool_obj[self._tops[mu]][1].jw_transform(self._fast_Pauli)
         Kmu_prev.mult_coeffs(self._pool_obj[self._tops[mu]][0])
 
         qc_psi.apply_operator(Kmu_prev)
@@ -166,7 +166,7 @@ def measure_gradient(self, params=None):
             else:
                 tamp = params[mu+1]
 
-            Kmu = self._pool_obj[self._tops[mu]][1].jw_transform()
+            Kmu = self._pool_obj[self._tops[mu]][1].jw_transform(self._fast_Pauli)
             Kmu.mult_coeffs(self._pool_obj[self._tops[mu]][0])
 
             Umu, pmu = trotterize(Kmu_prev, factor=-tamp, trotter_number=self._trotter_number)
@@ -212,7 +212,7 @@ def measure_gradient3(self):
         grads = np.zeros(len(self._pool_obj))
 
         for mu, (coeff, operator) in enumerate(self._pool_obj):
-            Kmu = operator.jw_transform()
+            Kmu = operator.jw_transform(self._fast_Pauli)
             Kmu.mult_coeffs(coeff)
             qc_psi.apply_operator(Kmu)
             grads[mu] = 2.0 * np.real(np.vdot(qc_sig.get_coeff_vec(), qc_psi.get_coeff_vec()))

src/qforte/adapters/molecule_adapters.py

@@ -38,6 +38,7 @@ def create_psi_mol(**kwargs):
     basis = kwargs['basis']
     multiplicity = kwargs['multiplicity']
     charge = kwargs['charge']
+    fast_Pauli = kwargs['fast_Pauli']
 
     qforte_mol = Molecule(mol_geometry = mol_geometry,
                                basis = basis,
@@ -164,6 +165,9 @@ def create_psi_mol(**kwargs):
                     mo_oeis[p, q] += 2 * mo_teis[p, q, i, i] - mo_teis[p, i, i, q]
 
     # Make hf_reference
+    #if fast_Pauli:
+    #    hf_reference = qforte.QubitBasis((1 << nel) - 1)
+    #else :
     hf_reference = [1] * (nel - 2 * frozen_core) + [0] * (2 * (nmo - frozen_virtual) - nel)
 
     # Build second quantized Hamiltonian
@@ -201,7 +205,7 @@ def create_psi_mol(**kwargs):
     qforte_mol.hf_energy = p4_Escf
     qforte_mol.hf_reference = hf_reference
     qforte_mol.sq_hamiltonian = Hsq
-    qforte_mol.hamiltonian = Hsq.jw_transform()
+    qforte_mol.hamiltonian = Hsq.jw_transform(fast_Pauli)
     qforte_mol.point_group = [point_group, irreps]
     qforte_mol.orb_irreps = orb_irreps
     qforte_mol.orb_irreps_to_int = orb_irreps_to_int
@@ -226,6 +230,8 @@ def create_external_mol(**kwargs):
         The qforte Molecule object which holds the molecular information.
     """
 
+    fast_Pauli = kwargs['fast_Pauli']
+
     qforte_mol = Molecule(multiplicity = kwargs['multiplicity'],
                                 charge = kwargs['charge'],
                                 filename = kwargs['filename'])
@@ -244,6 +250,9 @@ def create_external_mol(**kwargs):
     for p, q, r, s, h_pqrs in external_data['tei']['data']:
         qforte_sq_hamiltonian.add(h_pqrs/4.0, [p,q], [s,r]) # only works in C1 symmetry
 
+    #if fast_Pauli:
+    #    hf_reference = qforte.QubitBasis((1 << (external_data['na']['data'] + external_data['nb']['data'])) - 1)
+    #else:
     hf_reference = [0 for i in range(external_data['nso']['data'])]
     for n in range(external_data['na']['data'] + external_data['nb']['data']):
         hf_reference[n] = 1
@@ -252,6 +261,6 @@ def create_external_mol(**kwargs):
 
     qforte_mol.sq_hamiltonian = qforte_sq_hamiltonian
 
-    qforte_mol.hamiltonian = qforte_sq_hamiltonian.jw_transform()
+    qforte_mol.hamiltonian = qforte_sq_hamiltonian.jw_transform(fast_Pauli)
 
     return qforte_mol

src/qforte/bindings.cc

@@ -4,10 +4,13 @@
 
 #include "fmt/format.h"
 
+#include "helpers.h"
 #include "qubit_basis.h"
 #include "circuit.h"
 #include "gate.h"
 #include "computer.h"
+#include "pauli_string.h"
+#include "pauli_string_vector.h"
 #include "qubit_operator.h"
 #include "sq_operator.h"
 #include "sq_op_pool.h"
@@ -46,7 +49,8 @@ PYBIND11_MODULE(qforte, m) {
         .def("terms", &SQOperator::terms)
         .def("canonical_order", &SQOperator::canonical_order)
         .def("simplify", &SQOperator::simplify)
-        .def("jw_transform", &SQOperator::jw_transform)
+        .def("jw_transform", &SQOperator::jw_transform, py::arg("fast_Pauli") = true)
+        .def("fast_jw_transform", &SQOperator::fast_jw_transform)
         .def("str", &SQOperator::str)
         .def("__str__", &SQOperator::str)
         .def("__repr__", &SQOperator::str);
@@ -58,9 +62,9 @@ PYBIND11_MODULE(qforte, m) {
         .def("set_coeffs", &SQOpPool::set_coeffs)
         .def("terms", &SQOpPool::terms)
         .def("set_orb_spaces", &SQOpPool::set_orb_spaces)
-        .def("get_qubit_op_pool", &SQOpPool::get_qubit_op_pool)
+        .def("get_qubit_op_pool", &SQOpPool::get_qubit_op_pool, py::arg("fast_Pauli") = true)
         .def("get_qubit_operator", &SQOpPool::get_qubit_operator, py::arg("order_type"),
-             py::arg("combine_like_terms") = true)
+             py::arg("combine_like_terms") = true, py::arg("fast_Pauli") = true)
         .def("fill_pool", &SQOpPool::fill_pool)
         .def("str", &SQOpPool::str)
         .def("__getitem__", [](const SQOpPool &pool, size_t i) { return pool.terms()[i]; })
@@ -86,6 +90,7 @@ PYBIND11_MODULE(qforte, m) {
         .def("check_op_equivalence", &QubitOperator::check_op_equivalence)
         .def("num_qubits", &QubitOperator::num_qubits)
         .def("sparse_matrix", &QubitOperator::sparse_matrix)
+        .def("get_PauliStringVector", &QubitOperator::get_PauliStringVector)
         .def("str", &QubitOperator::str)
         .def("__str__", &QubitOperator::str)
         .def("__repr__", &QubitOperator::str);
@@ -111,13 +116,61 @@ PYBIND11_MODULE(qforte, m) {
     py::class_<QubitBasis>(m, "QubitBasis")
         .def(py::init<size_t>(), "n"_a = 0, "Make a basis element")
         .def("str", &QubitBasis::str)
-        .def("__str__", &QubitBasis::default_str)
-        .def("__repr__", &QubitBasis::default_str)
+        .def("__str__", [](const QubitBasis& qb) {
+            return qb.str(QubitBasis::max_qubits());
+        })
+        .def("__repr__", [](const QubitBasis& qb) {
+            return qb.str(QubitBasis::max_qubits());
+        })
         .def("flip_bit", &QubitBasis::flip_bit)
         .def("set_bit", &QubitBasis::set_bit)
-        .def("add", &QubitBasis::add)
+        .def("add", &QubitBasis::address)
+        .def("address", &QubitBasis::address)
         .def("get_bit", &QubitBasis::get_bit);
 
+    py::class_<PauliString>(m, "PauliString")
+        .def(py::init<size_t,size_t, std::complex<double>>(), "X"_a = 0,"Z"_a = 0, "coeff"_a = 1.0, "Make a Pauli string")
+        .def("str", &PauliString::str)
+        .def("__str__", &PauliString::str)
+        .def("__repr__", &PauliString::str)
+        .def("__eq__", [](const PauliString& lhs, const PauliString& rhs){return rhs == lhs;})
+        .def("__lt__", [](const PauliString& lhs, const PauliString& rhs){return lhs < rhs;})
+        .def("__mul__",[](const PauliString& lhs, const PauliString& rhs){return multiply(rhs,lhs);}, py::is_operator())
+        .def("__mul__",[](const PauliString& lhs, const std::complex<double> rhs){return multiply(lhs,rhs);}, py::is_operator())
+        .def("__rmul__",[](const PauliString& lhs, const std::complex<double> rhs){return multiply(lhs,rhs);}, py::is_operator())
+        .def("__add__",[](const PauliString& lhs, const PauliString& rhs){return add(lhs,rhs);}, py::is_operator())
+        .def("__sub__",[](const PauliString& lhs, const PauliString& rhs){return subtract(lhs,rhs);}, py::is_operator());
+
+   py::class_<PauliStringVector>(m, "PauliStringVector")
+        .def(py::init<std::vector<PauliString>>(), "Make a vector of Pauli string")
+        .def(py::init<>())
+        .def("str", &PauliStringVector::str)
+        .def("add_PauliString", &PauliStringVector::add_PauliString)
+        .def("add_PauliStringVector", &PauliStringVector::add_PauliStringVector)
+        .def("get_QubitOperator", &PauliStringVector::get_QubitOperator)
+        .def("num_qubits", &PauliStringVector::num_qubits)
+        .def("__str__", [](const PauliStringVector& vec){return join(vec.str(), "\n");})
+        .def("__repr__", [](const PauliStringVector& vec){return join(vec.str(), "\n");})
+        .def("__getitem__", [](const PauliStringVector& vec, unsigned int i){return vec[i];})
+        .def("__len__", [](const PauliStringVector& vec){return vec.get_vec().size();})
+        .def("__iter__", [](const PauliStringVector& PauliVector){
+                std::vector<PauliString> vec(PauliVector.get_vec());
+                return py::make_iterator(vec.begin(), vec.end());})
+        .def("__mul__",[](const PauliStringVector& lhs, const std::complex<double> rhs){return multiply(lhs,rhs);}, py::is_operator())
+        .def("__rmul__",[](const PauliStringVector& lhs, const std::complex<double> rhs){return multiply(lhs,rhs);}, py::is_operator())
+        .def("__mul__",[](const PauliStringVector& lhs, const PauliString& rhs){return multiply(lhs,rhs);}, py::is_operator())
+        .def("__rmul__",[](const PauliStringVector& lhs, const PauliString& rhs){return multiply(rhs,lhs);}, py::is_operator())
+        .def("__mul__",[](const PauliStringVector& lhs, const PauliStringVector& rhs){return multiply(rhs,lhs);}, py::is_operator())
+        .def("__add__",[](const PauliStringVector& lhs, const PauliString& rhs){return add(lhs,rhs);}, py::is_operator())
+        .def("__radd__",[](const PauliStringVector& lhs, const PauliString& rhs){return add(lhs,rhs);}, py::is_operator())
+        .def("__add__",[](const PauliStringVector& lhs, const PauliStringVector& rhs){return add(lhs,rhs);}, py::is_operator())
+        .def("__sub__",[](const PauliStringVector& lhs, const PauliString& rhs){return subtract(lhs,rhs);}, py::is_operator())
+        .def("__rsub__",[](const PauliStringVector& lhs, const PauliString& rhs){return subtract(rhs,lhs);}, py::is_operator())
+        .def("__sub__",[](const PauliStringVector& lhs, const PauliStringVector& rhs){return subtract(lhs,rhs);}, py::is_operator())
+        .def("get_vec", &PauliStringVector::get_vec)
+        .def("order_terms", &PauliStringVector::order_terms)
+        .def("simplify", &PauliStringVector::simplify);
+
     py::class_<Computer>(m, "Computer")
         .def(py::init<size_t>(), "nqubits"_a, "Make a quantum computer with 'nqubits' qubits")
         .def("apply_circuit_safe", &Computer::apply_circuit_safe)
@@ -185,6 +238,43 @@ PYBIND11_MODULE(qforte, m) {
         .def("reset", &local_timer::reset)
         .def("get", &local_timer::get);
 
+    m.def(
+        "exponentiate_fastpauli",
+        [](const PauliString& ps) {
+            if (std::abs(std::real(ps.coeff())) < 1.0e-14 and std::abs(std::imag(ps.coeff())) > 1.0e-14) {
+                PauliString ps1(0, 0, std::cos(std::imag(ps.coeff())));
+                PauliString ps2(ps.X(), ps.Z(), std::complex<double>{0,1} * std::sin(std::imag(ps.coeff())));
+                return add(ps1, ps2);}
+            std::string msg =
+                fmt::format("This function requires a pure imaginary coefficient. \n Invoked with ({:+f} {:+f} i)",
+                        std::real(ps.coeff()), std::imag(ps.coeff()));
+            throw std::invalid_argument(msg);
+        });
+
+    m.def(
+        "exponentiate_fastpauli",
+        [](const PauliStringVector& PauliVector) {
+            PauliString identity(0, 0, 1);
+            std::vector<PauliString> vec{identity};
+            PauliStringVector result(vec);
+            for(const PauliString& ps: PauliVector.get_vec()){
+                if (std::abs(std::real(ps.coeff())) > 1.0e-14 or std::abs(std::imag(ps.coeff())) < 1.0e-14) {
+                    std::string msg =
+                        fmt::format("This function requires a pure imaginary coefficient. \n Invoked with ({:+f} {:+f} i)",
+                            std::real(ps.coeff()), std::imag(ps.coeff()));
+                    throw std::invalid_argument(msg);
+                    }
+                else {
+                    PauliString ps1(0, 0, std::cos(std::imag(ps.coeff())));
+                    PauliString ps2(ps.X(), ps.Z(), std::complex<double>{0,1} * std::sin(std::imag(ps.coeff())));
+                    PauliStringVector ps_vector = add(ps1, ps2);
+                    result = multiply(ps_vector, result);
+                        }
+                    }
+            return result;
+                }
+        );
+
     m.def(
         "gate",
         [](std::string type, size_t target, std::complex<double> parameter) {

src/qforte/bitstring.cc

@@ -0,0 +1,42 @@
+#include "bitstring.h"
+
+std::string BitString::str(size_t nbit) const {
+    std::string s;
+    s += "|";
+    for (int i = 0; i < nbit; ++i) {
+        if (get_bit(i)) {
+            s += "1";
+        } else {
+            s += "0";
+        }
+    }
+    s += ">";
+    return s;
+}
+
+bool operator==(const BitString& lhs, const BitString& rhs){
+    return lhs.get_bits() == rhs.get_bits();
+}
+
+bool operator<(const BitString& lhs, const BitString& rhs){
+    return lhs.get_bits() < rhs.get_bits();
+}
+
+bool operator>(const BitString& lhs, const BitString& rhs){
+    return lhs.get_bits() > rhs.get_bits();
+}
+
+BitString operator^(const BitString& lhs, const BitString& rhs)
+{
+    return BitString(lhs.get_bits() ^ rhs.get_bits());
+}
+
+BitString operator&(const BitString& lhs, const BitString& rhs)
+{
+    return BitString(lhs.get_bits() & rhs.get_bits());
+}
+
+BitString operator|(const BitString& lhs, const BitString& rhs)
+{
+    return BitString(lhs.get_bits() | rhs.get_bits());
+}