Fermionic Excitation Implementation

src/tequila/quantumchemistry/chemistry_tools.py

@@ -2,12 +2,12 @@
 import typing
 import warnings
 from dataclasses import dataclass
-
+from copy import deepcopy
+from numbers import Real
 import numpy
 
-from tequila import BitString, QCircuit, TequilaException
+from tequila import BitString, QCircuit, TequilaException,Variable,compile_circuit
 from tequila.circuit import gates
-
 try:
     from openfermion.ops.representations import get_active_space_integrals  # needs openfermion 1.3
 except ImportError as E:
@@ -50,16 +50,132 @@ def __init__(self, generator, p0, transformation, indices=None, *args, **kwargs)
         self._name = "FermionicExcitation"
         self.transformation = transformation
         self.indices = indices
-
+        if isinstance(indices,tuple) and not hasattr(indices[0],"__len__"):
+            self.indices = [(indices[2 * i], indices[2 * i+1]) for i in range(len(indices) // 2)]
+        self.sign = self.format_excitation_variables(self.indices)
+        self.indices = self.format_excitation_indices(self.indices)
     def compile(self, *args, **kwargs):
         if self.is_convertable_to_qubit_excitation():
             target = []
             for x in self.indices:
                 for y in x:
                     target.append(y)
-            return gates.QubitExcitation(target=target, angle=-self.parameter, control=self.control)
+            return gates.QubitExcitation(target=target, angle=self.parameter, control=self.control)
+        else:
+            if self.transformation.lower().strip("_") == "jordanwigner":
+                return self.fermionic_excitation(angle=self.sign*self.parameter, indices=self.indices, control=self.control,opt=False)
+            else:
+                return gates.Trotterized(generator=self.generator, control=self.control, angle=self.parameter, steps=1)
+    def format_excitation_indices(self, idx):
+        """
+        Consistent formatting of excitation indices
+        idx = [(p0,q0),(p1,q1),...,(pn,qn)]
+        sorted as: p0<p1<pn and pi<qi
+        :param idx: list of index tuples describing a single(!) fermionic excitation
+        :return: list of index tuples
+        """
+
+        idx = [tuple(sorted(x)) for x in idx]
+        idx = sorted(idx, key=lambda x: x[0])
+        return list(idx)
+    def format_excitation_variables(self, idx):
+        """
+        Consistent formatting of excitation variable
+        idx = [(p0,q0),(p1,q1),...,(pn,qn)]
+        sorted as: pi<qi and p0 < p1 < p2
+        :param idx: list of index tuples describing a single(!) fermionic excitation
+        :return: sign of the variable with re-ordered indices
+        """
+        sig = 1
+        for pair in idx:
+            if pair[1]>pair[0]:
+                sig *= -1
+        for pair in range(len(idx)-1):
+            if idx[pair+1][0]>idx[pair][0]:
+                sig *= -1
+        return sig
+    def cCRy(self, target: int, dcontrol: typing.Union[list, int], control: typing.Union[list, int],
+             angle: typing.Union[Real, Variable, typing.Hashable], case: int = 1) -> QCircuit:
+        '''
+        Compilation of CRy as on https://doi.org/10.1103/PhysRevA.102.062612
+        If not control passed, Ry returned
+        Parameters
+        ----------
+        case: if 1 employs eq. 12 from the paper, if 0 eq. 13
+        '''
+        if control is not None and not len(control):
+            control = None
+        if isinstance(dcontrol, int):
+            dcontrol = [dcontrol]
+        if not len(dcontrol):
+            return compile_circuit(gates.Ry(angle=angle, target=target, control=control))
         else:
-            return gates.Trotterized(generator=self.generator, control=self.control, angle=self.parameter, steps=1)
+            if isinstance(angle, str):
+                angle = Variable(angle)
+            U = QCircuit()
+            aux = dcontrol[0]
+            ctr = deepcopy(dcontrol)
+            ctr.pop(0)
+            if case:
+                U += self.cCRy(target=target, dcontrol=ctr, angle=angle / 2, case=1, control=control) + gates.H(
+                    aux) + gates.CNOT(target, aux)
+                U += self.cCRy(target=target, dcontrol=ctr, angle=-angle / 2, case=0, control=control) + gates.CNOT(
+                    target, aux) + gates.H(aux)
+            else:
+                U += gates.H(aux) + gates.CNOT(target, aux) + self.cCRy(target=target, dcontrol=ctr, angle=-angle / 2,
+                                                                        case=0, control=control)
+                U += gates.CNOT(target, aux) + gates.H(aux) + self.cCRy(target=target, dcontrol=ctr, angle=angle / 2,
+                                                                        case=1, control=control)
+            return U
+
+    def fermionic_excitation(self, angle: typing.Union[Real, Variable, typing.Hashable], indices: typing.List,
+                             control: typing.Union[int, typing.List] = None, opt: bool = True) -> QCircuit:
+        '''
+            Excitation [(i,j),(k,l)],... compiled following https://doi.org/10.1103/PhysRevA.102.062612
+            opt: whether to optimized CNOT H CNOT --> Rz Rz CNOT Rz
+        '''
+        lto = []
+        lfrom = []
+        if isinstance(indices,tuple) and not hasattr(indices[0],"__len__"):
+            indices = [(indices[2 * i], indices[2 * i + 1]) for i in range(len(indices) // 2)]
+        for pair in indices:
+            lfrom.append(pair[0])
+            lto.append(pair[1])
+        Upair = QCircuit()
+        if isinstance(angle, str) or isinstance(angle, tuple):
+            angle = Variable(angle)
+        for i in range(len(lfrom) - 1):
+            Upair += gates.CNOT(lfrom[i + 1], lfrom[i])
+            Upair += gates.CNOT(lto[i + 1], lto[i])
+            Upair += gates.X(lto[i]) + gates.X(lfrom[i])
+        Upair += gates.CNOT(lto[-1], lfrom[-1])
+        crt = lfrom[::-1] + lto
+        Uladder = QCircuit()
+        pairs = lfrom + lto
+        pairs.sort()
+        orbs = []
+        for o in range(len(pairs) // 2):
+            orbs += [*range(pairs[2 * o] + 1, pairs[2 * o + 1])]
+        if len(orbs):
+            for o in range(len(orbs) - 1):
+                Uladder += gates.CNOT(orbs[o], orbs[o + 1])
+            Uladder += compile_circuit(gates.CZ(orbs[-1], lto[-1]))
+        crt.pop(-1)
+        if control is not None and (isinstance(control, int) or len(control) == 1):
+            if isinstance(control, int):
+                crt.append(control)
+            else:
+                crt = crt + control
+            control = []
+        Ur = self.cCRy(target=lto[-1], dcontrol=crt, angle=angle, control=control)
+        Upair2 = Upair.dagger()
+        if opt:
+            Ur.gates.pop(-1)
+            Ur.gates.pop(-1)
+            Upair2.gates.pop(0)
+            Ur += gates.Rz(numpy.pi / 2, target=lto[-1]) + gates.Rz(-numpy.pi / 2, target=lfrom[-1])
+            Ur += gates.CNOT(lto[-1], lfrom[-1]) + gates.Rz(numpy.pi / 2, target=lfrom[-1]) + gates.H(lfrom[-1])
+        return Upair + Uladder + Ur + Uladder.dagger() + Upair2
 
     def __str(self):
         if self.indices is not None:

src/tequila/quantumchemistry/qc_base.py

@@ -441,10 +441,14 @@ def make_excitation_gate(self, indices, angle, control=None, assume_real=True, *
 
         generator = self.make_excitation_generator(indices=indices, remove_constant_term=control is None)
         p0 = self.make_excitation_generator(indices=indices, form="P0", remove_constant_term=control is None)
-
+        if self.transformation.up_then_down:
+            idx = []
+            for pair in indices:
+                idx.append((pair[0]//2+(pair[0]%2)*self.n_orbitals,pair[1]//2+(pair[1]%2)*self.n_orbitals))
+        else:idx = indices
         return QCircuit.wrap_gate(
             FermionicGateImpl(angle=angle, generator=generator, p0=p0,
-                              transformation=type(self.transformation).__name__.lower(), indices=indices,
+                              transformation=type(self.transformation).__name__.lower(), indices=idx,
                               assume_real=assume_real,
                               control=control, **kwargs))
 

tests/test_chemistry.py

@@ -173,8 +173,6 @@ def test_ucc_singles_psi4():
 def do_test_ucc(qc_interface, parameters, result, trafo, backend="qulacs"):
     # check examples for comments
     psi4_interface = qc_interface(parameters=parameters, transformation=trafo)
-
-    hqc = psi4_interface.make_hamiltonian()
     amplitudes = psi4_interface.compute_ccsd_amplitudes()
     U = psi4_interface.make_uccsd_ansatz(trotter_steps=1, initial_amplitudes=amplitudes, include_reference_ansatz=True)
     variables = amplitudes.make_parameter_dictionary()
@@ -370,9 +368,9 @@ def test_hamiltonian_reduction(backend):
 
 @pytest.mark.skipif(condition=not HAS_PSI4 and not HAS_PYSCF, reason="psi4/pyscf not found")
 @pytest.mark.parametrize("assume_real", [True, False])
-@pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "tapered_bravyi_kitaev"])
+@pytest.mark.parametrize("trafo", ["jordan_wigner", "bravyi_kitaev", "reordered_jordan_wigner"])
 def test_fermionic_gates(assume_real, trafo):
-    mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.7\nLi 0.0 0.0 0.0", basis_set="sto-3g")
+    mol = tq.chemistry.Molecule(geometry="H 0.0 0.0 0.7\nLi 0.0 0.0 0.0", basis_set="sto-3g",transformation=trafo)
     U1 = mol.prepare_reference()
     U2 = mol.prepare_reference()
     variable_count = {}