Repository: quantumlib/OpenFermion
Pull request #54 (2017-10-17 00:07:45)

Move UCC independent functions back in OpenFermion
Moves the routines related to UCC that were independent of ProjectQ out of the OpenFermion-ProjectQ repo back into the main branch under utilities.  Merge this before merging branch.

babbush: You're good to merge when tests pass

On Mon, Oct 16, 2017 at 5:35 PM, Jarrod <notifications@github.com> wrote:

> *@jarrodmcc* commented on this pull request.
> ------------------------------
>
> In src/openfermion/utils/_unitary_cc.py
> <https://github.com/quantumlib/OpenFermion/pull/54#discussion_r145003841>:
>
> > +    spaces = range(n_virtual), range(n_occupied), range(2)
> +
> +    for i, j, s in itertools.product(*spaces):
> +        uccsd_generator += FermionOperator(
> +            (
> +                (2 * (i + n_occupied) + s, 1),
> +                (2 * j + s, 0),
> +            ),
> +            coefficient=t1[t1_ind(i, j)])
> +
> +        uccsd_generator += FermionOperator(
> +            (
> +                (2 * j + s, 1),
> +                (2 * (i + n_occupied) + s, 0),
> +            ),
> +            coefficient=-t1[t1_ind(i, j)])
>
> I removed the coefficient label since it wasn't needed. The reason it
> passed was that this is the PEP8 standard for keyword arguments (no space),
> even though it looks like an independent assignment line.
>
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> <https://github.com/quantumlib/OpenFermion/pull/54#discussion_r145003841>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/ANlTf57eYDn2V_IlXrRwzRVcYL8FA1Ktks5ss_Y7gaJpZM4P7Z55>
> .
>


Patch:
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index afc31ea..9b64c42 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -34,6 +34,11 @@ from ._sparse_tools import (expectation,
 
 from ._trotter_error import error_bound, error_operator
 
+from ._unitary_cc import (uccsd_convert_amplitude_format,
+                          uccsd_operator,
+                          uccsd_singlet_operator,
+                          uccsd_singlet_paramsize)
+
 # Imports out of alphabetical order to avoid circular dependancy.
 from ._dual_basis_trotter_error import (dual_basis_error_bound,
                                         dual_basis_error_operator)
diff --git a/src/openfermion/utils/_unitary_cc.py b/src/openfermion/utils/_unitary_cc.py
new file mode 100644
index 0000000..b0be28c
--- /dev/null
+++ b/src/openfermion/utils/_unitary_cc.py
@@ -0,0 +1,199 @@
+#   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.
+
+"""Module to create and manipulate unitary coupled cluster operators."""
+
+import itertools
+import numpy
+
+from openfermion.ops import FermionOperator, QubitOperator
+
+
+def uccsd_operator(single_amplitudes, double_amplitudes, anti_hermitian=True):
+    """Create a fermionic operator that is the generator of uccsd.
+
+    This a the most straight-forward method to generate UCCSD operators,
+    however it is slightly inefficient. In particular, it parameterizes
+    all possible excitations, so it represents a generalized unitary coupled
+    cluster ansatz, but also does not explicitly enforce the uniqueness
+    in parametrization, so it is redundant. For example there will be a linear
+    dependency in the ansatz of single_amplitudes[i,j] and
+    single_amplitudes[j,i].
+
+    Args:
+        single_amplitudes(list or ndarray): list of lists with each sublist
+            storing a list of indices followed by single excitation amplitudes
+            i.e. [[[i,j],t_ij], ...] OR [NxN] array storing single excitation
+            amplitudes corresponding to
+            t[i,j] * (a_i^\dagger a_j - H.C.)
+        double_amplitudes(list or ndarray): list of lists with each sublist
+            storing a list of indices followed by double excitation amplitudes
+            i.e. [[[i,j,k,l],t_ijkl], ...] OR [NxNxNxN] array storing double
+            excitation amplitudes corresponding to
+            t[i,j,k,l] * (a_i^\dagger a_j a_k^\dagger a_l - H.C.)
+        anti_hermitian(Bool): Flag to generate only normal CCSD operator
+            rather than unitary variant, primarily for testing
+
+    Returns:
+        uccsd_generator(FermionOperator): Anti-hermitian fermion operator that
+        is the generator for the uccsd wavefunction.
+    """
+
+    uccsd_generator = FermionOperator()
+
+    # Re-format inputs (ndarrays to lists) if necessary
+    if (isinstance(single_amplitudes, numpy.ndarray) or
+            isinstance(double_amplitudes, numpy.ndarray)):
+        single_amplitudes, double_amplitudes = uccsd_convert_amplitude_format(
+            single_amplitudes,
+            double_amplitudes)
+
+    # Add single excitations
+    for (i, j), t_ij in single_amplitudes:
+        i, j = int(i), int(j)
+        uccsd_generator += FermionOperator(((i, 1), (j, 0)), t_ij)
+        if anti_hermitian:
+            uccsd_generator += FermionOperator(((j, 1), (i, 0)), -t_ij)
+
+    # Add double excitations
+    for (i, j, k, l), t_ijkl in double_amplitudes:
+        i, j, k, l = int(i), int(j), int(k), int(l)
+        uccsd_generator += FermionOperator(
+            ((i, 1), (j, 0), (k, 1), (l, 0)), t_ijkl)
+        if anti_hermitian:
+            uccsd_generator += FermionOperator(
+                ((l, 1), (k, 0), (j, 1), (i, 0)), -t_ijkl)
+    return uccsd_generator
+
+
+def uccsd_convert_amplitude_format(single_amplitudes, double_amplitudes):
+    """Re-format single_amplitudes and double_amplitudes from ndarrays to lists.
+
+        Args:
+        single_amplitudes(ndarray): [NxN] array storing single excitation
+            amplitudes corresponding to t[i,j] * (a_i^\dagger a_j - H.C.)
+        double_amplitudes(ndarray): [NxNxNxN] array storing double excitation
+            amplitudes corresponding to
+            t[i,j,k,l] * (a_i^\dagger a_j a_k^\dagger a_l - H.C.)
+
+        Returns:
+        single_amplitudes_list(list): list of lists with each sublist storing
+            a list of indices followed by single excitation amplitudes
+            i.e. [[[i,j],t_ij], ...]
+        double_amplitudes_list(list): list of lists with each sublist storing
+            a list of indices followed by double excitation amplitudes
+            i.e. [[[i,j,k,l],t_ijkl], ...]
+    """
+    single_amplitudes_list, double_amplitudes_list = [], []
+
+    for i, j in zip(*single_amplitudes.nonzero()):
+        single_amplitudes_list.append([[i, j], single_amplitudes[i, j]])
+
+    for i, j, k, l in zip(*double_amplitudes.nonzero()):
+        double_amplitudes_list.append([[i, j, k, l],
+                                      double_amplitudes[i, j, k, l]])
+    return single_amplitudes_list, double_amplitudes_list
+
+
+def uccsd_singlet_paramsize(n_qubits, n_electrons):
+    """Determine number of independent amplitudes for singlet UCCSD
+
+    Args:
+        n_qubits(int): Number of qubits/spin-orbitals in the system
+        n_electrons(int): Number of electrons in the reference state
+
+    Returns:
+        Number of independent parameters for singlet UCCSD with a single
+        reference.
+    """
+    n_occupied = int(numpy.ceil(n_electrons / 2.))
+    n_virtual = n_qubits / 2 - n_occupied
+
+    n_single_amplitudes = n_occupied * n_virtual
+    n_double_amplitudes = n_single_amplitudes ** 2
+    return (n_single_amplitudes + n_double_amplitudes)
+
+
+def uccsd_singlet_operator(packed_amplitudes,
+                           n_qubits,
+                           n_electrons):
+    """Create a singlet UCCSD generator for a system with n_electrons
+
+    This function generates a FermionOperator for a UCCSD generator designed
+        to act on a single reference state consisting of n_qubits spin orbitals
+        and n_electrons electrons, that is a spin singlet operator, meaning it
+        conserves spin.
+
+    Args:
+        packed_amplitudes(ndarray): Compact array storing the unique single
+            and double excitation amplitudes for a singlet UCCSD operator.
+            The ordering lists unique single excitations before double
+            excitations.
+        n_qubits(int): Number of spin-orbitals used to represent the system,
+            which also corresponds to number of qubits in a non-compact map.
+        n_electrons(int): Number of electrons in the physical system.
+
+    Returns:
+        uccsd_generator(FermionOperator): Generator of the UCCSD operator that
+            builds the UCCSD wavefunction.
+    """
+    n_occupied = int(numpy.ceil(n_electrons / 2.))
+    n_virtual = int(n_qubits / 2 - n_occupied)  # Virtual Spatial Orbitals
+    n_t1 = int(n_occupied * n_virtual)
+
+    t1 = packed_amplitudes[:n_t1]
+    t2 = packed_amplitudes[n_t1:]
+
+    def t1_ind(i, j):
+        return i * n_occupied + j
+
+    def t2_ind(i, j, k, l):
+        return (i * n_occupied * n_virtual * n_occupied +
+                j * n_virtual * n_occupied +
+                k * n_occupied +
+                l)
+
+    uccsd_generator = FermionOperator()
+
+    # Define a compound space that is partitioned into occupied, virtual, spins
+    # the spin component assumes alpha spins are even, beta spins are odd
+    spaces = range(n_virtual), range(n_occupied), range(2)
+
+    # Generate all spin-conserving single excitations between occupied-virtual
+    for i, j, s in itertools.product(*spaces):
+        uccsd_generator += FermionOperator((
+            (2 * (i + n_occupied) + s, 1),
+            (2 * j + s, 0)),
+            t1[t1_ind(i, j)])
+
+        uccsd_generator += FermionOperator((
+            (2 * j + s, 1),
+            (2 * (i + n_occupied) + s, 0)),
+            -t1[t1_ind(i, j)])
+
+    # Generate all spin-conserving double excitations between occupied-virtual
+    for i, j, s, i2, j2, s2 in itertools.product(*spaces, repeat=2):
+        uccsd_generator += FermionOperator((
+            (2 * (i + n_occupied) + s, 1),
+            (2 * j + s, 0),
+            (2 * (i2 + n_occupied) + s2, 1),
+            (2 * j2 + s2, 0)),
+            t2[t2_ind(i, j, i2, j2)])
+
+        uccsd_generator += FermionOperator((
+            (2 * j2 + s2, 1),
+            (2 * (i2 + n_occupied) + s2, 0),
+            (2 * j + s, 1),
+            (2 * (i + n_occupied) + s, 0)),
+            -t2[t2_ind(i, j, i2, j2)])
+
+    return uccsd_generator


Test patch:
diff --git a/src/openfermion/utils/_unitary_cc_test.py b/src/openfermion/utils/_unitary_cc_test.py
new file mode 100644
index 0000000..e6f015c
--- /dev/null
+++ b/src/openfermion/utils/_unitary_cc_test.py
@@ -0,0 +1,205 @@
+#   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.
+
+"""Module to test unitary coupled cluster operators."""
+
+from __future__ import absolute_import
+
+import numpy
+import os
+import scipy
+import unittest
+
+from numpy.random import randn
+
+from openfermion.config import THIS_DIRECTORY
+from openfermion.hamiltonians import MolecularData
+from openfermion.ops import *
+from openfermion.transforms import *
+from openfermion.utils import *
+
+
+class UnitaryCC(unittest.TestCase):
+
+    def test_uccsd_anti_hermitian(self):
+        """Test operators are anti-Hermitian independent of inputs"""
+        test_orbitals = 4
+
+        single_amplitudes = randn(*(test_orbitals,) * 2)
+        double_amplitudes = randn(*(test_orbitals,) * 4)
+
+        generator = uccsd_operator(single_amplitudes, double_amplitudes)
+        conj_generator = hermitian_conjugated(generator)
+
+        self.assertTrue(generator.isclose(-1. * conj_generator))
+
+    def test_uccsd_singlet_anti_hermitian(self):
+        """Test that the singlet version is anti-Hermitian"""
+        test_orbitals = 8
+        test_electrons = 4
+
+        packed_amplitude_size = uccsd_singlet_paramsize(test_orbitals,
+                                                        test_electrons)
+
+        packed_amplitudes = randn(int(packed_amplitude_size))
+
+        generator = uccsd_singlet_operator(packed_amplitudes,
+                                           test_orbitals,
+                                           test_electrons)
+
+        conj_generator = hermitian_conjugated(generator)
+
+        self.assertTrue(generator.isclose(-1. * conj_generator))
+
+    def test_uccsd_singlet_build(self):
+        """Test a specific build of the UCCSD singlet operator"""
+        initial_amplitudes = [-1.14941450e-08, 5.65340614e-02]
+        n_orbitals = 4
+        n_electrons = 2
+
+        generator = uccsd_singlet_operator(initial_amplitudes,
+                                           n_orbitals,
+                                           n_electrons)
+
+        test_generator = (0.0565340614 * FermionOperator("2^ 0 3^ 1") +
+                          1.1494145e-08 * FermionOperator("1^ 3") +
+                          0.0565340614 * FermionOperator("3^ 1 2^ 0") +
+                          0.0565340614 * FermionOperator("2^ 0 2^ 0") +
+                          1.1494145e-08 * FermionOperator("0^ 2") +
+                          (-0.0565340614) * FermionOperator("1^ 3 0^ 2") +
+                          (-1.1494145e-08) * FermionOperator("3^ 1") +
+                          (-0.0565340614) * FermionOperator("1^ 3 1^ 3") +
+                          (-0.0565340614) * FermionOperator("0^ 2 0^ 2") +
+                          (-1.1494145e-08) * FermionOperator("2^ 0") +
+                          0.0565340614 * FermionOperator("3^ 1 3^ 1") +
+                          (-0.0565340614) * FermionOperator("0^ 2 1^ 3"))
+        self.assertTrue(test_generator.isclose(generator))
+
+    def test_sparse_uccsd_operator_numpy_inputs(self):
+        """Test numpy ndarray inputs to uccsd_operator that are sparse"""
+        test_orbitals = 30
+        sparse_single_amplitudes = numpy.zeros((test_orbitals, test_orbitals))
+        sparse_double_amplitudes = numpy.zeros((test_orbitals, test_orbitals,
+                                                test_orbitals, test_orbitals))
+
+        sparse_single_amplitudes[3, 5] = 0.12345
+        sparse_single_amplitudes[12, 4] = 0.44313
+
+        sparse_double_amplitudes[0, 12, 6, 2] = 0.3434
+        sparse_double_amplitudes[1, 4, 6, 13] = -0.23423
+
+        generator = uccsd_operator(sparse_single_amplitudes,
+                                   sparse_double_amplitudes)
+
+        test_generator = (0.12345 * FermionOperator("3^ 5") +
+                          (-0.12345) * FermionOperator("5^ 3") +
+                          0.44313 * FermionOperator("12^ 4") +
+                          (-0.44313) * FermionOperator("4^ 12") +
+                          0.3434 * FermionOperator("0^ 12 6^ 2") +
+                          (-0.3434) * FermionOperator("2^ 6 12^ 0") +
+                          (-0.23423) * FermionOperator("1^ 4 6^ 13") +
+                          0.23423 * FermionOperator("13^ 6 4^ 1"))
+        self.assertTrue(test_generator.isclose(generator))
+
+    def test_sparse_uccsd_operator_list_inputs(self):
+        """Test list inputs to uccsd_operator that are sparse"""
+        sparse_single_amplitudes = [[[3, 5], 0.12345],
+                                    [[12, 4], 0.44313]]
+        sparse_double_amplitudes = [[[0, 12, 6, 2], 0.3434],
+                                    [[1, 4, 6, 13], -0.23423]]
+
+        generator = uccsd_operator(sparse_single_amplitudes,
+                                   sparse_double_amplitudes)
+
+        test_generator = (0.12345 * FermionOperator("3^ 5") +
+                          (-0.12345) * FermionOperator("5^ 3") +
+                          0.44313 * FermionOperator("12^ 4") +
+                          (-0.44313) * FermionOperator("4^ 12") +
+                          0.3434 * FermionOperator("0^ 12 6^ 2") +
+                          (-0.3434) * FermionOperator("2^ 6 12^ 0") +
+                          (-0.23423) * FermionOperator("1^ 4 6^ 13") +
+                          0.23423 * FermionOperator("13^ 6 4^ 1"))
+        self.assertTrue(test_generator.isclose(generator))
+
+    def test_ucc(self):
+        geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
+        basis = 'sto-3g'
+        multiplicity = 1
+        filename = os.path.join(THIS_DIRECTORY, 'data',
+                                'H2_sto-3g_singlet_0.7414')
+        self.molecule = MolecularData(
+            geometry, basis, multiplicity, filename=filename)
+        self.molecule.load()
+
+        # Get molecular Hamiltonian.
+        self.molecular_hamiltonian = self.molecule.get_molecular_hamiltonian()
+
+        # Get FCI RDM.
+        self.fci_rdm = self.molecule.get_molecular_rdm(use_fci=1)
+
+        # Get explicit coefficients.
+        self.nuclear_repulsion = self.molecular_hamiltonian.constant
+        self.one_body = self.molecular_hamiltonian.one_body_tensor
+        self.two_body = self.molecular_hamiltonian.two_body_tensor
+
+        # Get fermion Hamiltonian.
+        self.fermion_hamiltonian = normal_ordered(get_fermion_operator(
+            self.molecular_hamiltonian))
+
+        # Get qubit Hamiltonian.
+        self.qubit_hamiltonian = jordan_wigner(self.fermion_hamiltonian)
+
+        # Get the sparse matrix.
+        self.hamiltonian_matrix = get_sparse_operator(
+            self.molecular_hamiltonian)
+        # Test UCCSD for accuracy against FCI using loaded t amplitudes.
+        ucc_operator = uccsd_operator(
+            self.molecule.ccsd_single_amps,
+            self.molecule.ccsd_double_amps)
+
+        hf_state = jw_hartree_fock_state(
+            self.molecule.n_electrons, count_qubits(self.qubit_hamiltonian))
+        uccsd_sparse = jordan_wigner_sparse(ucc_operator)
+        uccsd_state = scipy.sparse.linalg.expm_multiply(uccsd_sparse,
+                                                        hf_state)
+        expected_uccsd_energy = expectation(self.hamiltonian_matrix,
+                                            uccsd_state)
+        self.assertAlmostEqual(expected_uccsd_energy, self.molecule.fci_energy,
+                               places=4)
+        print("UCCSD ENERGY: {}".format(expected_uccsd_energy))
+
+        # Test CCSD for precise match against FCI using loaded t amplitudes.
+        ccsd_operator = uccsd_operator(
+            self.molecule.ccsd_single_amps,
+            self.molecule.ccsd_double_amps,
+            anti_hermitian=False)
+
+        ccsd_sparse_r = jordan_wigner_sparse(ccsd_operator)
+        ccsd_sparse_l = jordan_wigner_sparse(
+            -hermitian_conjugated(ccsd_operator))
+
+        # Test CCSD for precise match against FCI using loaded t amplitudes
+        ccsd_operator = uccsd_operator(
+            self.molecule.ccsd_single_amps,
+            self.molecule.ccsd_double_amps,
+            anti_hermitian=False)
+
+        ccsd_sparse_r = jordan_wigner_sparse(ccsd_operator)
+        ccsd_sparse_l = jordan_wigner_sparse(
+            -hermitian_conjugated(ccsd_operator))
+        ccsd_state_r = scipy.sparse.linalg.expm_multiply(ccsd_sparse_r,
+                                                         hf_state)
+        ccsd_state_l = scipy.sparse.linalg.expm_multiply(ccsd_sparse_l,
+                                                         hf_state)
+        expected_ccsd_energy = ccsd_state_l.getH().dot(
+            self.hamiltonian_matrix.dot(ccsd_state_r))[0, 0]
+        self.assertAlmostEqual(expected_ccsd_energy, self.molecule.fci_energy)