Reduced BCS type Hamiltonians (quantumlib#770)

* improved import formatting

* support RichardsonGaudin-type Hamiltonians

* improved formatting

* corrected docstring which corresponds to the ladder operator

Co-authored-by: Nicholas Rubin <[email protected]>

src/openfermion/_compat_test.py

@@ -16,12 +16,12 @@
 import pytest
 import deprecation
 
-import openfermion
-from openfermion._compat import wrap_module
-
 from cirq._compat import deprecated
 from cirq.testing import assert_deprecated
 
+import openfermion
+from openfermion._compat import wrap_module
+
 
 def deprecated_test(test: Callable) -> Callable:
     """Marks a test as using deprecated functionality.

src/openfermion/hamiltonians/__init__.py

@@ -51,6 +51,8 @@
 
 from .mean_field_dwave import mean_field_dwave
 
+from .richardson_gaudin import RichardsonGaudin
+
 from .plane_wave_hamiltonian import (
     dual_basis_external_potential,
     plane_wave_external_potential,

src/openfermion/hamiltonians/richardson_gaudin.py

@@ -0,0 +1,109 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""This module constructs Hamiltonians of the Richardson Gaudin type.
+"""
+from itertools import chain, product
+import numpy
+from openfermion.ops.representations import (PolynomialTensor,
+                                             get_tensors_from_integrals)
+from openfermion.ops.representations import DOCIHamiltonian
+from openfermion.ops import QubitOperator
+
+
+class RichardsonGaudin(DOCIHamiltonian):
+    r"""Richardson Gaudin model.
+
+    Class for storing and constructing Richardson Gaudin hamiltonians
+    combining an equi-distant potential ladder like potential per
+    qubit with a uniform coupling between any pair of
+    qubits with coupling strength g, which can be either attractive
+    (g<0) or repulsive (g>0).
+
+    The operators represented by this class has the form:
+
+    .. math::
+
+        H = \sum_{p=0} (p + 1) N_p + g/2 \sum_{p < q} P_p^\dagger P_q,
+
+    where
+
+    .. math::
+
+        \begin{align}
+        N_p &= (1 - \sigma^Z_p)/2, \\
+        P_p &= a_{p,\beta} a_{p,\alpha} = S^{-} = \sigma^X + i \sigma^Y, \\
+        g &= constant coupling term
+        \end{align}
+
+    Note;
+        The diagonal of the Hamiltonian is composed of the values in
+        range((n_qubits+1)*n_qubits//2+1).
+    """
+
+    def __init__(self, g, n_qubits):
+        r"""Richardson Gaudin model on a given number of qubits.
+
+        Args:
+            g (float): Coupling strength
+            n_qubits (int): Number of qubits
+        """
+        hc = numpy.zeros((n_qubits,))
+        hr1 = numpy.zeros((n_qubits, n_qubits))
+        hr2 = numpy.zeros((n_qubits, n_qubits))
+        for p in range(n_qubits):
+            hc[p] = 2 * (p + 1)
+            for q in range(n_qubits):
+                if p != q:
+                    hr1[p, q] = g
+        super().__init__(0, hc, hr1, hr2)
+
+    @DOCIHamiltonian.constant.setter
+    def constant(self, value):
+        raise TypeError('Raw edits of the constant of a RichardsonGaudin model'
+                        'is not allowed. Either adjust the g parameter '
+                        'or cast to another PolynomialTensor class.')
+
+    @DOCIHamiltonian.n_body_tensors.setter
+    def n_body_tensors(self, value):
+        raise TypeError(
+            'Raw edits of the n_body_tensors of a RichardsonGaudin model'
+            'is not allowed. Either adjust the g parameter '
+            'or cast to another PolynomialTensor class.')
+
+    def get_antisymmetrized_tensors(self):
+        r"""Antisymmetrized Tensors
+        Directly returns antisymmetrized tensors, which, when used
+        to construct an FermionOperator via an InteractionOperator
+        produce a FermionOperator that acts like this RichardsonGaudin
+        Hamiltonian on the paired (seniority zero) subspace.
+        Compared to the FermionOperator that can be obtained via the
+        n_body_tensors property from the DOCIHamiltonian class
+        the FermionOperator from the tensors returned by this function
+        do not contain same spin coupling terms. These terms
+        act trivially on the paired subspace and this the two Hamiltonian
+        agree on any senioirty zero state.
+        Returns:
+            tuple: Tuple of one and two body tensors.
+        """
+        g = self.hr1[0, 1]
+        spatial_orbs = self.hc.shape[0]
+        h1 = numpy.diag(numpy.arange(spatial_orbs) + 1)
+        h1 = numpy.kron(h1, numpy.eye(2))
+        h2 = numpy.zeros((2 * spatial_orbs,) * 4)
+        for p, q in product(range(spatial_orbs), repeat=2):
+            if p != q:
+                h2[2 * p, 2 * p + 1, 2 * q + 1, 2 * q] = g / 2
+                h2[2 * p + 1, 2 * p, 2 * q, 2 * q + 1] = g / 2
+
+        h2 = h2 - numpy.einsum('ijlk', h2)
+
+        return h1, h2

src/openfermion/hamiltonians/richardson_gaudin_test.py

@@ -0,0 +1,92 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""Tests for Richardson Gaudin model module."""
+
+import pytest
+import numpy as np
+
+from openfermion.hamiltonians import RichardsonGaudin
+from openfermion.ops import QubitOperator
+from openfermion.transforms import get_fermion_operator
+from openfermion.transforms import jordan_wigner
+from openfermion.linalg import get_sparse_operator
+from openfermion import(get_fermion_operator, InteractionOperator, \
+                         normal_ordered)
+
+
+@pytest.mark.parametrize('g, n_qubits, expected', [
+    (0.3, 2,
+     QubitOperator('3.0 [] + 0.15 [X0 X1] + \
+0.15 [Y0 Y1] - 1.0 [Z0] - 2.0 [Z1]')),
+    (-0.1, 3,
+     QubitOperator('6.0 [] - 0.05 [X0 X1] - 0.05 [X0 X2] - \
+0.05 [Y0 Y1] - 0.05 [Y0 Y2] - 1.0 [Z0] - 0.05 [X1 X2] - \
+0.05 [Y1 Y2] - 2.0 [Z1] - 3.0 [Z2]')),
+])
+def test_richardson_gaudin_hamiltonian(g, n_qubits, expected):
+    rg = RichardsonGaudin(g, n_qubits)
+    rg_qubit = rg.qubit_operator
+    assert rg_qubit == expected
+
+    assert np.array_equal(
+        np.sort(np.unique(get_sparse_operator(rg_qubit).diagonal())),
+        2 * np.array(list(range((n_qubits + 1) * n_qubits // 2 + 1))))
+
+
+def test_n_body_tensor_errors():
+    rg = RichardsonGaudin(1.7, n_qubits=2)
+    with pytest.raises(TypeError):
+        rg.n_body_tensors = 0
+    with pytest.raises(TypeError):
+        rg.constant = 1.1
+
+
+@pytest.mark.parametrize("g, n_qubits", [(0.2, 4), (-0.2, 4)])
+def test_fermionic_hamiltonian_from_integrals(g, n_qubits):
+    rg = RichardsonGaudin(g, n_qubits)
+    #hc, hr1, hr2 = rg.hc, rg.hr1, rg.hr2
+    doci = rg
+    constant = doci.constant
+    reference_constant = 0
+
+    doci_qubit_op = doci.qubit_operator
+    doci_mat = get_sparse_operator(doci_qubit_op).toarray()
+    doci_eigvals = np.linalg.eigh(doci_mat)[0]
+
+    tensors = doci.n_body_tensors
+    one_body_tensors, two_body_tensors = tensors[(1, 0)], tensors[(1, 1, 0, 0)]
+
+    fermion_op = get_fermion_operator(
+        InteractionOperator(constant, one_body_tensors, 0.5 * two_body_tensors))
+    fermion_op = normal_ordered(fermion_op)
+    fermion_mat = get_sparse_operator(fermion_op).toarray()
+    fermion_eigvals = np.linalg.eigh(fermion_mat)[0]
+
+    one_body_tensors2, two_body_tensors2 = rg.get_antisymmetrized_tensors()
+    fermion_op2 = get_fermion_operator(
+        InteractionOperator(reference_constant, one_body_tensors2,
+                            0.5 * two_body_tensors2))
+    fermion_op2 = normal_ordered(fermion_op2)
+    fermion_mat2 = get_sparse_operator(fermion_op2).toarray()
+    fermion_eigvals2 = np.linalg.eigh(fermion_mat2)[0]
+
+    for eigval in doci_eigvals:
+        assert any(abs(fermion_eigvals -
+                       eigval) < 1e-6), "The DOCI spectrum should have \
+        been contained in the spectrum of the fermionic operator constructed via the \
+        DOCIHamiltonian class"
+
+    for eigval in doci_eigvals:
+        assert any(abs(fermion_eigvals2 -
+                       eigval) < 1e-6), "The DOCI spectrum should have \
+       been contained in the spectrum of the fermionic operators constructed via the anti \
+       symmetrized tensors"

src/openfermion/ops/operators/symbolic_operator.py

@@ -121,7 +121,9 @@ def different_indices_commute(self):
 
     def __init__(self, term=None, coefficient=1.):
         if not isinstance(coefficient, COEFFICIENT_TYPES):
-            raise ValueError('Coefficient must be a numeric type.')
+            raise ValueError(
+                'Coefficient must be a numeric type. Got {}'.format(
+                    type(coefficient)))
 
         # Initialize the terms dictionary
         self.terms = {}

src/openfermion/ops/representations/doci_hamiltonian.py

@@ -125,7 +125,7 @@ def z_term(self, p):
             [QubitOperator] -- Z term on the chosen qubit.
         """
         return QubitOperator("Z" + str(p),
-                             self.hc[p] / 2 + sum(self.hr2[:, p]) / 2)
+                             -self.hc[p] / 2 - sum(self.hr2[:, p]) / 2)
 
     def zz_term(self, p, q):
         """Returns the ZZ term on a single pair of qubits as a QubitOperator
@@ -357,9 +357,53 @@ def zero(cls, n_qubits):
                    numpy.zeros((n_qubits,) * 2, dtype=numpy.complex128),
                    numpy.zeros((n_qubits,) * 2, dtype=numpy.complex128))
 
+    def get_projected_integrals(self):
+        ''' Creates the one and two body integrals that would correspond to a
+        hypothetic electronic structure Hamiltonian, which would satisfy the
+        given set of hc, hr1 and hr2.
+
+        This is technically not well-defined, as hr2 is
+        not generated in a one-to-one fashion. This implies that calling
+
+        get_doci_from_integrals(
+            *get_projected_integrals_from_doci(
+               hc, hr1, hr2
+            )
+        )
+
+        should return the same hc, hr1, and hr2, but there is no such guarantee
+        for
+
+        get_projected_integrals_from_doci(
+           *get_doci_from_integrals(
+              one_body_integrals, two_body_integrals
+           )
+        )
+
+        but this method attemps to create integrals that conform to the
+        same symmetries as a physical electronic structure Hamiltonian would,
+        with inevitable loss of information due to the ambiguity above.
+
+        Args:
+           hc [numpy array]: The single-particle DOCI terms in matrix form
+           hr1 [numpy array]: The off-diagonal DOCI Hamiltonian terms in matrix
+               form
+           hr2 [numpy array]: The diagonal DOCI Hamiltonian terms in matrix form
+
+        Returns:
+           projected_onebody_integrals [numpy array]: The corresponding one-body
+               integrals for the electronic structure Hamiltonian
+           projected_twobody_integrals [numpy array]: The corresponding two body
+               integrals for the electronic structure Hamiltonian
+        '''
+        one_body_integrals, two_body_integrals = \
+            get_projected_integrals_from_doci(self.hc, self.hr1, self.hr2)
+        return one_body_integrals, two_body_integrals
+
 
 def get_tensors_from_doci(hc, hr1, hr2):
-    '''Makes the one and two-body tensors from the DOCI Hamiltonian matrices
+    '''Makes the one and two-body tensors of a fermionic "parent Hamoltonian" \
+    from the DOCI Hamiltonian matrices
 
     Args:
         hc [numpy array]: The single-particle DOCI terms in matrix form
@@ -377,6 +421,10 @@ def get_tensors_from_doci(hc, hr1, hr2):
         get_projected_integrals_from_doci(hc, hr1, hr2)
     one_body_coefficients, two_body_coefficients = get_tensors_from_integrals(
         one_body_integrals, two_body_integrals)
+
+    two_body_coefficients = two_body_coefficients - numpy.einsum(
+        'ijlk', two_body_coefficients)
+
     return one_body_coefficients, two_body_coefficients
 
 
@@ -410,17 +458,26 @@ def get_projected_integrals_from_doci(hc, hr1, hr2):
             integrals for the electronic structure Hamiltonian
     '''
     n_qubits = hr1.shape[0]
-    projected_onebody_integrals = numpy.zeros((n_qubits, n_qubits))
+    projected_onebody_integrals = numpy.zeros((n_qubits, n_qubits),
+                                              dtype=hc.dtype)
     projected_twobody_integrals = numpy.zeros(
-        (n_qubits, n_qubits, n_qubits, n_qubits))
+        (n_qubits, n_qubits, n_qubits, n_qubits), dtype=hc.dtype)
     for p in range(n_qubits):
         projected_onebody_integrals[p, p] = hc[p] / 2
         projected_twobody_integrals[p, p, p, p] = hr2[p, p]
         for q in range(n_qubits):
-            if p == q:
+            if p <= q:
                 continue
-            projected_twobody_integrals[p, q, q, p] = hr2[p, q] / 2
-            projected_twobody_integrals[p, p, q, q] = hr1[p, q]
+
+            projected_twobody_integrals[p, q, q,
+                                        p] = hr2[p, q] / 2 + hr1[p, q] / 2
+            projected_twobody_integrals[q, p, p,
+                                        q] = hr2[q, p] / 2 + hr1[p, q] / 2
+
+            projected_twobody_integrals[p, p, q, q] += hr1[p, q]
+            projected_twobody_integrals[p, q, p, q] += hr1[p, q]
+            projected_twobody_integrals[q, q, p, p] += hr1[p, q]
+            projected_twobody_integrals[q, p, q, p] += hr1[p, q]
 
     return projected_onebody_integrals, projected_twobody_integrals