pUCCD ansatz

tangelo/algorithms/variational/vqe_solver.py

@@ -46,6 +46,7 @@ class BuiltInAnsatze(Enum):
     VSQS = agen.VSQS
     UCCGD = agen.UCCGD
     ILC = agen.ILC
+    pUCCD = agen.pUCCD
 
 
 class VQESolver:
@@ -125,6 +126,9 @@ def __init__(self, opt_dict):
                 warnings.warn("Spin-orbital ordering shifted to all spin-up first then down to ensure efficient generator screening "
                               "for the Jordan-Wigner mapping with QCC-based ansatze.", RuntimeWarning)
                 self.up_then_down = True
+            if self.ansatz == BuiltInAnsatze.pUCCD and self.qubit_mapping.lower() != "hcb":
+                warnings.warn("Forcing the hard-core boson mapping for the pUCCD ansatz.", RuntimeWarning)
+                self.mapping = "HCB"
             # QCC and QMF and ILC require a reference state that can be represented by a single layer of RZ-RX gates on each qubit.
             # This decomposition can not be determined from a general Circuit reference state.
             if isinstance(self.ref_state, Circuit):
@@ -210,6 +214,8 @@ def build(self):
             if isinstance(self.ansatz, BuiltInAnsatze):
                 if self.ansatz in {BuiltInAnsatze.UCC1, BuiltInAnsatze.UCC3}:
                     self.ansatz = self.ansatz.value
+                elif self.ansatz == BuiltInAnsatze.pUCCD:
+                    self.ansatz = self.ansatz.value(self.molecule, **self.ansatz_options)
                 elif self.ansatz in self.builtin_ansatze:
                     self.ansatz = self.ansatz.value(self.molecule, self.qubit_mapping, self.up_then_down, **self.ansatz_options)
                 else:

tangelo/toolboxes/ansatz_generator/__init__.py

@@ -24,3 +24,4 @@
 from .hea import HEA
 from .variational_circuit import VariationalCircuitAnsatz
 from .uccgd import UCCGD
+from .puccd import pUCCD

tangelo/toolboxes/ansatz_generator/ansatz_utils.py

@@ -425,10 +425,10 @@ def givens_gate(target, theta, is_variational=False):
     """Generates the list of gates corresponding to a givens rotation exp(-theta*(XX+YY))
 
     Explicitly the two-qubit matrix is
-    [[0,      0,           0,       0],
+    [[1,      0,           0,       0],
      [0,  cos(theta), -sin(theta),  0],
      [0,  sin(theta),  cos(theta),  0],
-     [0,      0,            0,      0]]
+     [0,      0,            0,      1]]
 
     Args:
         target (list): list of two integers that indicate which qubits are involved in the givens rotation

tangelo/toolboxes/ansatz_generator/puccd.py

@@ -0,0 +1,205 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# 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 defines the pUCCD ansatz class. The molecular FermionOperator is
+expected to be converted to a BosonOperator (electrons in pairs). Single bosonic
+excitations (corresponding to double fermion excitations) form the ansatz. Those
+excitations are transformed into a quantum circuit via Givens rotations.
+
+Ref:
+    - V.E. Elfving, M. Millaruelo, J.A. Gámez, and C. Gogolin.
+        Phys. Rev. A 103, 032605 (2021).
+"""
+
+import itertools
+import numpy as np
+
+from tangelo.linq import Circuit
+from tangelo.toolboxes.ansatz_generator import Ansatz
+from tangelo.toolboxes.ansatz_generator.ansatz_utils import givens_gate
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import vector_to_circuit
+
+
+class pUCCD(Ansatz):
+    """This class implements the pUCCD ansatz, as described in Phys. Rev. A 103,
+    032605 (2021). Electrons are described as hard-core boson and only double
+    excitations are considered.
+
+    Args:
+        molecule (SecondQuantizedMolecule): Self-explanatory.
+        reference_state (string): String refering to an initial state.
+            Default: "HF".
+    """
+
+    def __init__(self, molecule, reference_state="HF"):
+
+        if molecule.spin != 0:
+            raise NotImplementedError("pUCCD is implemented only for closed-shell system.")
+
+        self.molecule = molecule
+        self.n_spatial_orbitals = molecule.n_active_mos
+        self.n_electrons = molecule.n_active_electrons
+
+        # Set total number of parameters.
+        self.n_occupied = int(self.n_electrons / 2)
+        self.n_virtual = self.n_spatial_orbitals - self.n_occupied
+        self.n_var_params = self.n_occupied * self.n_virtual
+
+        # Supported reference state initialization.
+        self.supported_reference_state = {"HF", "zero"}
+        # Supported var param initialization
+        self.supported_initial_var_params = {"zeros", "ones", "random"}
+
+        # Default initial parameters for initialization.
+        self.var_params_default = "ones"
+        self.reference_state = reference_state
+
+        self.var_params = None
+        self.circuit = None
+
+    def set_var_params(self, var_params=None):
+        """Set values for variational parameters, such as ones, zeros or random
+        numbers providing some keywords for users, and also supporting direct
+        user input (list or numpy array). Return the parameters so that
+        workflows such as VQE can retrieve these values.
+        """
+        if var_params is None:
+            var_params = self.var_params_default
+
+        if isinstance(var_params, str):
+            var_params = var_params.lower()
+            if (var_params not in self.supported_initial_var_params):
+                raise ValueError(f"Supported keywords for initializing variational parameters: {self.supported_initial_var_params}")
+            if var_params == "ones":
+                initial_var_params = np.ones((self.n_var_params,), dtype=float)
+            elif var_params == "random":
+                initial_var_params = 2.e-1 * (np.random.random((self.n_var_params,)) - 0.5)
+        else:
+            initial_var_params = np.array(var_params)
+            if initial_var_params.size != self.n_var_params:
+                raise ValueError(f"Expected {self.n_var_params} variational parameters but "
+                                 f"received {initial_var_params.size}.")
+        self.var_params = initial_var_params
+        return initial_var_params
+
+    def prepare_reference_state(self):
+        """Returns circuit preparing the reference state of the ansatz (e.g
+        prepare reference wavefunction with HF, multi-reference state, etc).
+        These preparations must be consistent with the transform used to obtain
+        the qubit operator.
+        """
+
+        if self.reference_state not in self.supported_reference_state:
+            raise ValueError(f"Only supported reference state methods are:{self.supported_reference_state}")
+
+        if self.reference_state == "HF":
+            vector = [1 if i < self.n_electrons // 2 else 0 for i in range(self.n_spatial_orbitals)]
+            return vector_to_circuit(vector)
+        elif self.reference_state == "zero":
+            return Circuit()
+
+    def build_circuit(self, var_params=None):
+        """Build and return the quantum circuit implementing the state
+        preparation ansatz (with currently specified initial_state and
+        var_params).
+        """
+
+        if var_params is not None:
+            self.set_var_params(var_params)
+        elif self.var_params is None:
+            self.set_var_params()
+
+        excitations = self._get_double_excitations()
+
+        # Prepend reference state circuit
+        reference_state_circuit = self.prepare_reference_state()
+
+        # Obtain quantum circuit through trivial trotterization of the qubit operator
+        # Keep track of the order in which pauli words have been visited for fast subsequent parameter updates
+        self.exc_to_param_mapping = dict()
+        rotation_gates = []
+
+        # Parallel ordering (rotations on different qubits can happen at the
+        # same time.
+        excitations_per_layer = [[]]
+        free_qubits_per_layer = [set(range(self.n_spatial_orbitals))]
+
+        # Classify excitations into circuit layers (single pass on all
+        # excitations).
+        for p, q in excitations:
+            excitations_added = False
+            for qubit_indices, gates in zip(free_qubits_per_layer, excitations_per_layer):
+                if p in qubit_indices and q in qubit_indices:
+                    qubit_indices -= {p, q}
+                    gates += [(p, q)]
+                    excitations_added = True
+                    break
+
+            # If the excitation cannot be added to at least one previous layer,
+            # create a new layer.
+            if not excitations_added:
+                excitations_per_layer.append([(p, q)])
+                qubit_indices = set(range(self.n_spatial_orbitals))
+                qubit_indices -= {p, q}
+                free_qubits_per_layer.append(qubit_indices)
+
+        excitations = list(itertools.chain.from_iterable(excitations_per_layer))
+        self.exc_to_param_mapping = {v: k for k, v in enumerate(excitations)}
+
+        rotation_gates = [givens_gate((p, q), 0., is_variational=True) for (p, q) in excitations]
+        rotation_gates = list(itertools.chain.from_iterable(rotation_gates))
+
+        puccd_circuit = Circuit(rotation_gates)
+
+        # Skip over the reference state circuit if it is empty.
+        if reference_state_circuit.size != 0:
+            self.circuit = reference_state_circuit + puccd_circuit
+        else:
+            self.circuit = puccd_circuit
+
+        self.update_var_params(self.var_params)
+        return self.circuit
+
+    def update_var_params(self, var_params):
+        """Shortcut: set value of variational parameters in the already-built
+        ansatz circuit member. Preferable to rebuilt your circuit from scratch,
+        which can be an involved process.
+        """
+
+        self.set_var_params(var_params)
+        var_params = self.var_params
+
+        excitations = self._get_double_excitations()
+        for i, (p, q) in enumerate(excitations):
+            gate_index = self.exc_to_param_mapping[(p, q)]
+            self.circuit._variational_gates[gate_index].parameter = var_params[i]
+
+    def _get_double_excitations(self):
+        """Construct the UCC double excitations for the given amount of occupied
+        and virtual orbitals.
+
+        Returns:
+            list of int tuples: List of (p, q) excitations corresponding to the
+                occupied orbital p to virtual orbital q.
+        """
+
+        # Generate double indices in seniority 0 space.
+        excitations = list()
+        for p, q in itertools.product(
+            range(self.n_occupied),
+            range(self.n_occupied, self.n_occupied+self.n_virtual)
+        ):
+            excitations += [(p, q)]
+
+        return excitations

tangelo/toolboxes/ansatz_generator/tests/test_puccd.py

@@ -0,0 +1,73 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# 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.
+
+import unittest
+import numpy as np
+
+from tangelo.molecule_library import mol_H2_sto3g
+from tangelo.toolboxes.qubit_mappings.hcb import hard_core_boson_operator, boson_to_qubit_mapping
+from tangelo.toolboxes.ansatz_generator.puccd import pUCCD
+from tangelo.linq import get_backend
+
+
+class pUCCDTest(unittest.TestCase):
+
+    def test_puccd_set_var_params(self):
+        """Verify behavior of set_var_params for different inputs (keyword,
+        list, numpy array).
+        """
+
+        puccd_ansatz = pUCCD(mol_H2_sto3g)
+
+        one_ones = np.ones((1,))
+
+        puccd_ansatz.set_var_params("ones")
+        np.testing.assert_array_almost_equal(puccd_ansatz.var_params, one_ones, decimal=6)
+
+        puccd_ansatz.set_var_params([1.])
+        np.testing.assert_array_almost_equal(puccd_ansatz.var_params, one_ones, decimal=6)
+
+        puccd_ansatz.set_var_params(np.array([1.]))
+        np.testing.assert_array_almost_equal(puccd_ansatz.var_params, one_ones, decimal=6)
+
+    def test_puccd_incorrect_number_var_params(self):
+        """Return an error if user provide incorrect number of variational
+        parameters.
+        """
+
+        puccd_ansatz = pUCCD(mol_H2_sto3g)
+
+        self.assertRaises(ValueError, puccd_ansatz.set_var_params, np.array([1., 1., 1., 1.]))
+
+    def test_puccd_H2(self):
+        """Verify closed-shell pUCCD functionalities for H2."""
+
+        # Build circuit.
+        puccd_ansatz = pUCCD(mol_H2_sto3g)
+        puccd_ansatz.build_circuit()
+
+        # Build qubit hamiltonian for energy evaluation.
+        qubit_hamiltonian = boson_to_qubit_mapping(
+            hard_core_boson_operator(mol_H2_sto3g.fermionic_hamiltonian)
+        )
+
+        # Assert energy returned is as expected for given parameters.
+        sim = get_backend()
+        puccd_ansatz.update_var_params([-0.22617753])
+        energy = sim.get_expectation_value(qubit_hamiltonian, puccd_ansatz.circuit)
+        self.assertAlmostEqual(energy, -1.13727, delta=1e-4)
+
+
+if __name__ == "__main__":
+    unittest.main()

tangelo/toolboxes/molecular_computation/coefficients.py

@@ -0,0 +1,51 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# 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 containing functions to manipulate molecular coefficient arrays."""
+
+import numpy as np
+
+
+def spatial_from_spinorb(one_body_coefficients, two_body_coefficients):
+    """Function to reverse openfermion.chem.molecular_data.spinorb_from_spatial.
+
+    Args:
+        one_body_coefficients: One-body coefficients (array of 2N*2N, where N
+            is the number of molecular orbitals).
+        two_body_coefficients: Two-body coefficients (array of 2N*2N*2N*2N,
+            where N is the number of molecular orbitals).
+
+    Returns:
+        (array of floats, array of floats): One- and two-body integrals (arrays
+            of N*N and N*N*N*N elements, where N is the number of molecular
+            orbitals.
+    """
+    # Get the number of MOs = number of SOs / 2.
+    n_mos = one_body_coefficients.shape[0] // 2
+
+    # Initialize Hamiltonian integrals.
+    one_body_integrals = np.zeros((n_mos, n_mos), dtype=complex)
+    two_body_integrals = np.zeros((n_mos, n_mos, n_mos, n_mos), dtype=complex)
+
+    # Loop through coefficients.
+    for p in range(n_mos):
+        for q in range(n_mos):
+            # Populate 1-body integrals.
+            one_body_integrals[p, q] = one_body_coefficients[2*p, 2*q]
+            # Continue looping to prepare 2-body integrals.
+            for r in range(n_mos):
+                for s in range(n_mos):
+                    two_body_integrals[p, q, r, s] = two_body_coefficients[2*p, 2*q+1, 2*r+1, 2*s]
+
+    return one_body_integrals, two_body_integrals