Repository: quantumlib/OpenFermion
Pull request #522 (2019-05-22 00:58:11)

JW for complex coefficients


Patch:
diff --git a/src/openfermion/ops/_interaction_operator.py b/src/openfermion/ops/_interaction_operator.py
index d1be1aa..6349aca 100644
--- a/src/openfermion/ops/_interaction_operator.py
+++ b/src/openfermion/ops/_interaction_operator.py
@@ -121,7 +121,8 @@ class InteractionOperator(PolynomialTensor):
 
     @classmethod
     def zero(cls, n_qubits):
-        return cls(0, numpy.zeros((n_qubits,) * 2), numpy.zeros((n_qubits,) * 4))
+        return cls(0, numpy.zeros((n_qubits,) * 2, dtype=numpy.complex128),
+                   numpy.zeros((n_qubits,) * 4, dtype=numpy.complex128))
 
     def projected(self, indices, exact=False):
         projected_n_body_tensors = self.projected_n_body_tensors(
diff --git a/src/openfermion/transforms/_jordan_wigner.py b/src/openfermion/transforms/_jordan_wigner.py
index 9592f6b..44095f5 100644
--- a/src/openfermion/transforms/_jordan_wigner.py
+++ b/src/openfermion/transforms/_jordan_wigner.py
@@ -152,7 +152,7 @@ def jordan_wigner_interaction_op(iop, n_qubits=None):
     # Transform other one-body terms and "diagonal" two-body terms
     for p, q in itertools.combinations(range(n_qubits), 2):
         # One-body
-        coefficient = .5 * (iop[(p, 1), (q, 0)] + iop[(q, 1), (p, 0)])
+        coefficient = .5 * (iop[(p, 1), (q, 0)] + iop[(q, 1), (p, 0)].conjugate())
         qubit_operator += jordan_wigner_one_body(p, q, coefficient)
 
         # Two-body
@@ -166,17 +166,22 @@ def jordan_wigner_interaction_op(iop, n_qubits=None):
     for (p, q), (r, s) in itertools.combinations(
             itertools.combinations(range(n_qubits), 2),
             2):
-        coefficient = (iop[(p, 1), (q, 1), (r, 0), (s, 0)] -
-                       iop[(p, 1), (q, 1), (s, 0), (r, 0)] -
-                       iop[(q, 1), (p, 1), (r, 0), (s, 0)] +
-                       iop[(q, 1), (p, 1), (s, 0), (r, 0)])
+        coefficient = 0.5 * (
+                iop[(p, 1), (q, 1), (r, 0), (s, 0)] +
+                iop[(s, 1), (r, 1), (q, 0), (p, 0)].conjugate() -
+                iop[(p, 1), (q, 1), (s, 0), (r, 0)] -
+                iop[(r, 1), (s, 1), (q, 0), (p, 0)].conjugate() -
+                iop[(q, 1), (p, 1), (r, 0), (s, 0)] -
+                iop[(s, 1), (r, 1), (p, 0), (q, 0)].conjugate() +
+                iop[(q, 1), (p, 1), (s, 0), (r, 0)] +
+                iop[(r, 1), (s, 1), (p, 0), (q, 0)].conjugate())
         qubit_operator += jordan_wigner_two_body(p, q, r, s, coefficient)
 
     return qubit_operator
 
 
 def jordan_wigner_one_body(p, q, coefficient=1.):
-    """Map the term a^\dagger_p a_q + a^\dagger_q a_p to QubitOperator.
+    """Map the term a^\dagger_p a_q + h.c. to QubitOperator.
 
     Note that the diagonal terms are divided by a factor of 2
     because they are equal to their own Hermitian conjugate.
@@ -184,11 +189,16 @@ def jordan_wigner_one_body(p, q, coefficient=1.):
     # Handle off-diagonal terms.
     qubit_operator = QubitOperator()
     if p != q:
-        a, b = sorted([p, q])
-        parity_string = tuple((z, 'Z') for z in range(a + 1, b))
-        for operator in ['X', 'Y']:
-            operators = ((a, operator),) + parity_string + ((b, operator),)
-            qubit_operator += QubitOperator(operators, .5 * coefficient)
+        if p > q:
+            p, q = q, p
+            coefficient = coefficient.conjugate()
+        parity_string = tuple((z, 'Z') for z in range(p + 1, q))
+        for c, (op_a, op_b) in [(coefficient.real, 'XX'),
+                                (coefficient.real, 'YY'),
+                                (coefficient.imag, 'YX'),
+                                (-coefficient.imag, 'XY')]:
+            operators = ((p, op_a),) + parity_string + ((q, op_b),)
+            qubit_operator += QubitOperator(operators, .5 * c)
 
     # Handle diagonal terms.
     else:
@@ -213,23 +223,28 @@ def jordan_wigner_two_body(p, q, r, s, coefficient=1.):
 
     # Handle case of four unique indices.
     elif len(set([p, q, r, s])) == 4:
-
+        if (p > q) ^ (r > s):
+            coefficient *= -1
         # Loop through different operators which act on each tensor factor.
-        for operator_p, operator_q, operator_r in itertools.product(
-                ['X', 'Y'], repeat=3):
-            if [operator_p, operator_q, operator_r].count('X') % 2:
-                operator_s = 'X'
+        for ops in itertools.product('XY', repeat=4):
+            op_p, op_q, op_r, op_s = ops
+            # Get coefficients.
+            if ops.count('X') % 2:
+                coeff = .125 * coefficient.imag
+                if ''.join(ops) in ['XYXX', 'YXXX', 'YYXY', 'YYYX']:
+                    coeff *= -1
             else:
-                operator_s = 'Y'
+                coeff = .125 * coefficient.real
+                if ''.join(ops) not in ['XXYY', 'YYXX']:
+                    coeff *= -1
+            if not coeff:
+                continue
 
             # Sort operators.
             [(a, operator_a), (b, operator_b),
-             (c, operator_c), (d, operator_d)] = sorted(
-                [(p, operator_p), (q, operator_q),
-                 (r, operator_r), (s, operator_s)],
-                key=lambda pair: pair[0])
+             (c, operator_c), (d, operator_d)] = sorted(zip([p, q, r, s], ops))
 
-            # Computer operator strings.
+            # Compute operator strings.
             operators = ((a, operator_a),)
             operators += tuple((z, 'Z') for z in range(a + 1, b))
             operators += ((b, operator_b),)
@@ -237,16 +252,6 @@ def jordan_wigner_two_body(p, q, r, s, coefficient=1.):
             operators += tuple((z, 'Z') for z in range(c + 1, d))
             operators += ((d, operator_d),)
 
-            # Get coefficients.
-            coeff = .125 * coefficient
-            parity_condition = bool(operator_p != operator_q or
-                                    operator_p == operator_r)
-            if (p > q) ^ (r > s):
-                if not parity_condition:
-                    coeff *= -1.
-            elif parity_condition:
-                coeff *= -1.
-
             # Add term.
             qubit_operator += QubitOperator(operators, coeff)
 
@@ -255,32 +260,49 @@ def jordan_wigner_two_body(p, q, r, s, coefficient=1.):
 
         # Identify equal tensor factors.
         if p == r:
-            a, b = sorted([q, s])
+            if q > s:
+                a, b = s, q
+                coefficient = -coefficient.conjugate()
+            else:
+                a, b = q, s
+                coefficient = -coefficient
             c = p
         elif p == s:
-            a, b = sorted([q, r])
+            if q > r:
+                a, b = r, q
+                coefficient = coefficient.conjugate()
+            else:
+                a, b = q, r
             c = p
         elif q == r:
-            a, b = sorted([p, s])
+            if p > s:
+                a, b = s, p
+                coefficient = coefficient.conjugate()
+            else:
+                a, b = p, s
             c = q
         elif q == s:
-            a, b = sorted([p, r])
+            if p > r:
+                a, b = r, p
+                coefficient = -coefficient.conjugate()
+            else:
+                a, b = p, r
+                coefficient = -coefficient
             c = q
 
         # Get operators.
         parity_string = tuple((z, 'Z') for z in range(a + 1, b))
         pauli_z = QubitOperator(((c, 'Z'),))
-        for operator in ['X', 'Y']:
-            operators = ((a, operator),) + parity_string + ((b, operator),)
-
-            # Get coefficient.
-            if (p == s) or (q == r):
-                coeff = .25 * coefficient
-            else:
-                coeff = -.25 * coefficient
+        for c, (op_a, op_b) in [(coefficient.real, 'XX'),
+                                (coefficient.real, 'YY'),
+                                (coefficient.imag, 'YX'),
+                                (-coefficient.imag, 'XY')]:
+            operators = ((a, op_a),) + parity_string + ((b, op_b),)
+            if not c:
+                continue
 
             # Add term.
-            hopping_term = QubitOperator(operators, coeff)
+            hopping_term = QubitOperator(operators, c / 4)
             qubit_operator -= pauli_z * hopping_term
             qubit_operator += hopping_term
 


Test patch:
diff --git a/src/openfermion/transforms/_jordan_wigner_test.py b/src/openfermion/transforms/_jordan_wigner_test.py
index 977c58e..705765c 100644
--- a/src/openfermion/transforms/_jordan_wigner_test.py
+++ b/src/openfermion/transforms/_jordan_wigner_test.py
@@ -12,6 +12,7 @@
 
 """Tests  _jordan_wigner.py."""
 from __future__ import absolute_import
+import itertools
 import os
 import unittest
 
@@ -226,7 +227,7 @@ class InteractionOperatorsJWTest(unittest.TestCase):
         """Test consistency with JW for FermionOperators."""
         # Random interaction operator
         n_qubits = 5
-        iop = random_interaction_operator(n_qubits)
+        iop = random_interaction_operator(n_qubits, real=False)
         op1 = jordan_wigner(iop)
         op2 = jordan_wigner(get_fermion_operator(iop))
 
@@ -252,42 +253,43 @@ class InteractionOperatorsJWTest(unittest.TestCase):
         # Make sure it agrees with jordan_wigner(FermionTerm).
         for p in range(self.n_qubits):
             for q in range(self.n_qubits):
+                coefficient = numpy.random.randn()
+                if p != q:
+                    coefficient += 1j * numpy.random.randn()
                 # Get test qubit operator.
-                test_operator = jordan_wigner_one_body(p, q)
+                test_operator = jordan_wigner_one_body(p, q, coefficient)
 
                 # Get correct qubit operator.
-                fermion_term = FermionOperator(((p, 1), (q, 0)))
+                fermion_term = FermionOperator(((p, 1), (q, 0)), coefficient)
+                if p != q:
+                    fermion_term += FermionOperator(((q, 1), (p, 0)),
+                            coefficient.conjugate())
                 correct_op = jordan_wigner(fermion_term)
-                hermitian_conjugate = hermitian_conjugated(fermion_term)
-                if not fermion_term == hermitian_conjugate:
-                    correct_op += jordan_wigner(hermitian_conjugate)
 
                 self.assertTrue(test_operator == correct_op)
 
     def test_jordan_wigner_two_body(self):
         # Make sure it agrees with jordan_wigner(FermionTerm).
-        for p in range(self.n_qubits):
-            for q in range(self.n_qubits):
-                for r in range(self.n_qubits):
-                    for s in range(self.n_qubits):
-                        # Get test qubit operator.
-                        test_operator = jordan_wigner_two_body(p, q, r, s)
-
-                        # Get correct qubit operator.
-                        fermion_term = FermionOperator(((p, 1), (q, 1),
-                                                        (r, 0), (s, 0)))
-                        correct_op = jordan_wigner(fermion_term)
-                        hermitian_conjugate = hermitian_conjugated(
-                            fermion_term)
-                        if not fermion_term == hermitian_conjugate:
-                            if p == r and q == s:
-                                pass
-                            else:
-                                correct_op += jordan_wigner(
-                                    hermitian_conjugate)
-
-                        self.assertTrue(test_operator == correct_op,
-                                        str(test_operator - correct_op))
+        for p, q, r, s in itertools.product(range(self.n_qubits), repeat=4):
+            coefficient = numpy.random.randn()
+            if set([p, q]) != set([r, s]):
+                coefficient += 1j * numpy.random.randn()
+
+            # Get test qubit operator.
+            test_operator = jordan_wigner_two_body(
+                    p, q, r, s, coefficient)
+
+            # Get correct qubit operator.
+            fermion_term = FermionOperator(((p, 1), (q, 1),
+                                            (r, 0), (s, 0)), coefficient)
+            if set([p, q]) != set([r, s]):
+                fermion_term += FermionOperator(
+                        ((s, 1), (r, 1), (q, 0), (p, 0)),
+                        coefficient.conjugate())
+            correct_op = jordan_wigner(fermion_term)
+
+            self.assertTrue(test_operator == correct_op,
+                            str(test_operator - correct_op))
 
     def test_jordan_wigner_twobody_interaction_op_allunique(self):
         test_op = FermionOperator('1^ 2^ 3 4')