# Source code for pennylane.ops.qubit.qchem_ops

# Copyright 2018-2021 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""
This submodule contains the discrete-variable quantum operations that come
from quantum chemistry applications.
"""
# pylint:disable=abstract-method,arguments-differ,protected-access
import cmath
import math
import numpy as np

import pennylane as qml
from pennylane.operation import Operation

INV_SQRT2 = 1 / math.sqrt(2)

# Four term gradient recipe for controlled rotations
c1 = INV_SQRT2 * (np.sqrt(2) + 1) / 4
c2 = INV_SQRT2 * (np.sqrt(2) - 1) / 4
a = np.pi / 2
b = 3 * np.pi / 2
four_term_grad_recipe = ([[c1, 1, a], [-c1, 1, -a], [-c2, 1, b], [c2, 1, -b]],)

[docs]class SingleExcitation(Operation):
r"""SingleExcitation(phi, wires)
Single excitation rotation.

.. math:: U(\phi) = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & \cos(\phi/2) & -\sin(\phi/2) & 0 \\
0 & \sin(\phi/2) & \cos(\phi/2) & 0 \\
0 & 0 & 0 & 1
\end{bmatrix}.

This operation performs a rotation in the two-dimensional subspace :math:\{|01\rangle,
|10\rangle\}. The name originates from the occupation-number representation of
fermionic wavefunctions, where the transformation  from :math:|10\rangle to :math:|01\rangle
is interpreted as "exciting" a particle from the first qubit to the second.

**Details:**

* Number of wires: 2
* Number of parameters: 1
* Gradient recipe: The SingleExcitation operator satisfies a four-term parameter-shift rule
(see Appendix F, https://arxiv.org/abs/2104.05695)

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int]): the wires the operation acts on

**Example**

The following circuit performs the transformation :math:|10\rangle\rightarrow \cos(
\phi/2)|10\rangle -\sin(\phi/2)|01\rangle:

.. code-block::

dev = qml.device('default.qubit', wires=2)

@qml.qnode(dev)
def circuit(phi):
qml.PauliX(wires=0)
qml.SingleExcitation(phi, wires=[0, 1])
return qml.state()

circuit(0.1)
"""

num_params = 1
num_wires = 2
par_domain = "R"
generator = [
np.array([[0, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, 0]]),
-1 / 2,
]
has_unitary_generator = False

@classmethod
def _matrix(cls, *params):
theta = params
c = math.cos(theta / 2)
s = math.sin(theta / 2)

return np.array([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1]])

[docs]    @staticmethod
def decomposition(theta, wires):
decomp_ops = [
qml.CNOT(wires=[wires, wires]),
qml.CRY(theta, wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
]
return decomp_ops

(phi,) = self.parameters
return SingleExcitation(-phi, wires=self.wires)

[docs]class SingleExcitationMinus(Operation):
r"""SingleExcitationMinus(phi, wires)
Single excitation rotation with negative phase-shift outside the rotation subspace.

.. math:: U_-(\phi) = \begin{bmatrix}
e^{-i\phi/2} & 0 & 0 & 0 \\
0 & \cos(\phi/2) & -\sin(\phi/2) & 0 \\
0 & \sin(\phi/2) & \cos(\phi/2) & 0 \\
0 & 0 & 0 & e^{-i\phi/2}
\end{bmatrix}.

**Details:**

* Number of wires: 2
* Number of parameters: 1
* Gradient recipe: :math:\frac{d}{d\phi}f(U_-(\phi)) = \frac{1}{2}\left[f(U_-(\phi+\pi/2)) - f(U_-(\phi-\pi/2))\right]
where :math:f is an expectation value depending on :math:U_-(\phi).

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int] or int): the wires the operation acts on

"""
num_params = 1
num_wires = 2
par_domain = "R"
generator = [
np.array([[1, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]]),
-1 / 2,
]
has_unitary_generator = True

@classmethod
def _matrix(cls, *params):
theta = params
c = math.cos(theta / 2)
s = math.sin(theta / 2)
e = cmath.exp(-1j * theta / 2)

return np.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]])

[docs]    @staticmethod
def decomposition(theta, wires):
decomp_ops = [
qml.PauliX(wires=wires),
qml.PauliX(wires=wires),
qml.ControlledPhaseShift(-theta / 2, wires=[wires, wires]),
qml.PauliX(wires=wires),
qml.PauliX(wires=wires),
qml.ControlledPhaseShift(-theta / 2, wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.CRY(theta, wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
]
return decomp_ops

(phi,) = self.parameters
return SingleExcitationMinus(-phi, wires=self.wires)

[docs]class SingleExcitationPlus(Operation):
r"""SingleExcitationPlus(phi, wires)
Single excitation rotation with positive phase-shift outside the rotation subspace.

.. math:: U_+(\phi) = \begin{bmatrix}
e^{i\phi/2} & 0 & 0 & 0 \\
0 & \cos(\phi/2) & -\sin(\phi/2) & 0 \\
0 & \sin(\phi/2) & \cos(\phi/2) & 0 \\
0 & 0 & 0 & e^{i\phi/2}
\end{bmatrix}.

**Details:**

* Number of wires: 2
* Number of parameters: 1
* Gradient recipe: :math:\frac{d}{d\phi}f(U_+(\phi)) = \frac{1}{2}\left[f(U_+(\phi+\pi/2)) - f(U_+(\phi-\pi/2))\right]
where :math:f is an expectation value depending on :math:U_+(\phi).

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int] or int): the wires the operation acts on

"""
num_params = 1
num_wires = 2
par_domain = "R"
generator = [
np.array([[-1, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, -1]]),
-1 / 2,
]
has_unitary_generator = True

@classmethod
def _matrix(cls, *params):
theta = params
c = math.cos(theta / 2)
s = math.sin(theta / 2)
e = cmath.exp(1j * theta / 2)

return np.array([[e, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, e]])

[docs]    @staticmethod
def decomposition(theta, wires):
decomp_ops = [
qml.PauliX(wires=wires),
qml.PauliX(wires=wires),
qml.ControlledPhaseShift(theta / 2, wires=[wires, wires]),
qml.PauliX(wires=wires),
qml.PauliX(wires=wires),
qml.ControlledPhaseShift(theta / 2, wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.CRY(theta, wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
]
return decomp_ops

(phi,) = self.parameters
return SingleExcitationPlus(-phi, wires=self.wires)

[docs]class DoubleExcitation(Operation):
r"""DoubleExcitation(phi, wires)
Double excitation rotation.

This operation performs an :math:SO(2) rotation in the two-dimensional subspace :math:\{
|1100\rangle,|0011\rangle\}. More precisely, it performs the transformation

.. math::

&|0011\rangle \rightarrow \cos(\phi/2) |0011\rangle - \sin(\phi/2) |1100\rangle\\
&|1100\rangle \rightarrow \cos(\phi/2) |1100\rangle + \sin(\phi/2) |0011\rangle,

while leaving all other basis states unchanged.

The name originates from the occupation-number representation of fermionic wavefunctions, where
the transformation from :math:|1100\rangle to :math:|0011\rangle is interpreted as
"exciting" two particles from the first pair of qubits to the second pair of qubits.

**Details:**

* Number of wires: 4
* Number of parameters: 1
* Gradient recipe: The DoubleExcitation operator satisfies a four-term parameter-shift rule
(see Appendix F, https://arxiv.org/abs/2104.05695):

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int]): the wires the operation acts on

**Example**

The following circuit performs the transformation :math:|1100\rangle\rightarrow \cos(
\phi/2)|1100\rangle +\sin(\phi/2)|0011\rangle):

.. code-block::

dev = qml.device('default.qubit', wires=4)

@qml.qnode(dev)
def circuit(phi):
qml.PauliX(wires=0)
qml.PauliX(wires=1)
qml.DoubleExcitation(phi, wires=[0, 1, 2, 3])
return qml.state()

circuit(0.1)
"""

num_params = 1
num_wires = 4
par_domain = "R"

G = np.zeros((16, 16), dtype=np.complex64)
G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
generator = [G, -1 / 2]
has_unitary_generator = False

@classmethod
def _matrix(cls, *params):
theta = params
c = math.cos(theta / 2)
s = math.sin(theta / 2)

U = np.eye(16)
U[3, 3] = c  # 3 (dec) = 0011 (bin)
U[3, 12] = -s  # 12 (dec) = 1100 (bin)
U[12, 3] = s
U[12, 12] = c

return U

[docs]    @staticmethod
def decomposition(theta, wires):
# This decomposition is the "upside down" version of that on p17 of https://arxiv.org/abs/2104.05695
decomp_ops = [
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.RY(theta / 8, wires=wires),
qml.RY(-theta / 8, wires=wires),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.RY(theta / 8, wires=wires),
qml.RY(-theta / 8, wires=wires),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.RY(-theta / 8, wires=wires),
qml.RY(theta / 8, wires=wires),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.RY(-theta / 8, wires=wires),
qml.RY(theta / 8, wires=wires),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
]

return decomp_ops

(theta,) = self.parameters
return DoubleExcitation(-theta, wires=self.wires)

[docs]class DoubleExcitationPlus(Operation):
r"""DoubleExcitationPlus(phi, wires)
Double excitation rotation with positive phase-shift outside the rotation subspace.

This operation performs an :math:SO(2) rotation in the two-dimensional subspace :math:\{
|1100\rangle,|0011\rangle\} while applying a phase-shift on other states. More precisely,
it performs the transformation

.. math::

&|0011\rangle \rightarrow \cos(\phi/2) |0011\rangle - \sin(\phi/2) |1100\rangle\\
&|1100\rangle \rightarrow \cos(\phi/2) |1100\rangle + \sin(\phi/2) |0011\rangle\\
&|x\rangle \rightarrow e^{i\phi/2} |x\rangle,

for all other basis states :math:|x\rangle.

**Details:**

* Number of wires: 4
* Number of parameters: 1
* Gradient recipe: :math:\frac{d}{d\phi}f(U_+(\phi)) = \frac{1}{2}\left[f(U_+(\phi+\pi/2)) - f(U_+(\phi-\pi/2))\right]
where :math:f is an expectation value depending on :math:U_+(\phi)

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int]): the wires the operation acts on
"""

num_params = 1
num_wires = 4
par_domain = "R"

G = -1 * np.eye(16, dtype=np.complex64)
G[3, 3] = 0
G[12, 12] = 0
G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
generator = [G, -1 / 2]
has_unitary_generator = True

@classmethod
def _matrix(cls, *params):
theta = params
c = math.cos(theta / 2)
s = math.sin(theta / 2)
e = cmath.exp(1j * theta / 2)

U = e * np.eye(16, dtype=np.complex64)
U[3, 3] = c  # 3 (dec) = 0011 (bin)
U[3, 12] = -s  # 12 (dec) = 1100 (bin)
U[12, 3] = s
U[12, 12] = c

return U

(theta,) = self.parameters
return DoubleExcitationPlus(-theta, wires=self.wires)

[docs]class DoubleExcitationMinus(Operation):
r"""DoubleExcitationMinus(phi, wires)
Double excitation rotation with negative phase-shift outside the rotation subspace.

This operation performs an :math:SO(2) rotation in the two-dimensional subspace :math:\{
|1100\rangle,|0011\rangle\} while applying a phase-shift on other states. More precisely,
it performs the transformation

.. math::

&|0011\rangle \rightarrow \cos(\phi/2) |0011\rangle - \sin(\phi/2) |1100\rangle\\
&|1100\rangle \rightarrow \cos(\phi/2) |1100\rangle + \sin(\phi/2) |0011\rangle\\
&|x\rangle \rightarrow e^{-i\phi/2} |x\rangle,

for all other basis states :math:|x\rangle.

**Details:**

* Number of wires: 4
* Number of parameters: 1
* Gradient recipe: :math:\frac{d}{d\phi}f(U_-(\phi)) = \frac{1}{2}\left[f(U_-(\phi+\pi/2)) - f(U_-(\phi-\pi/2))\right]
where :math:f is an expectation value depending on :math:U_-(\phi)

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int]): the wires the operation acts on
"""

num_params = 1
num_wires = 4
par_domain = "R"

G = np.eye(16, dtype=np.complex64)
G[3, 3] = 0
G[12, 12] = 0
G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
G[12, 3] = 1j  # 12 (dec) = 1100 (bin)
generator = [G, -1 / 2]
has_unitary_generator = True

@classmethod
def _matrix(cls, *params):
theta = params
c = math.cos(theta / 2)
s = math.sin(theta / 2)
e = cmath.exp(-1j * theta / 2)

U = e * np.eye(16, dtype=np.complex64)
U[3, 3] = c  # 3 (dec) = 0011 (bin)
U[3, 12] = -s  # 12 (dec) = 1100 (bin)
U[12, 3] = s
U[12, 12] = c

return U

(theta,) = self.parameters
return DoubleExcitationMinus(-theta, wires=self.wires)

[docs]class OrbitalRotation(Operation):
r"""OrbitalRotation(phi, wires)

For two neighbouring spatial orbitals :math:\{|\Phi_{0}\rangle, |\Phi_{1}\rangle\}, this operation
performs the following transformation

.. math::
&|\Phi_{0}\rangle = \cos(\phi/2)|\Phi_{0}\rangle - \sin(\phi/2)|\Phi_{1}\rangle\\
&|\Phi_{1}\rangle = \cos(\phi/2)|\Phi_{0}\rangle + \sin(\phi/2)|\Phi_{1}\rangle,

with the same orbital operation applied in the :math:\alpha and :math:\beta spin orbitals.

.. figure:: ../../_static/qchem/orbital_rotation_decomposition_extended.png
:align: center
:width: 100%
:target: javascript:void(0);

Here, :math:G(\phi) represents a single-excitation Givens rotation, implemented in PennyLane
as the :class:~.SingleExcitation operation.

**Details:**

* Number of wires: 4
* Number of parameters: 1
* Gradient recipe: The OrbitalRotation operator satisfies the four-term parameter-shift rule
(see Appendix F, https://arxiv.org/abs/2104.05695)

Args:
phi (float): rotation angle :math:\phi
wires (Sequence[int]): the wires the operation acts on

**Example**

.. code-block::

>>> dev = qml.device('default.qubit', wires=4)
>>> @qml.qnode(dev)
... def circuit(phi):
...     qml.BasisState(np.array([1, 1, 0, 0]), wires=[0, 1, 2, 3])
...     qml.OrbitalRotation(phi, wires=[0, 1, 2, 3])
...     return qml.state()
>>> circuit(0.1)
array([ 0.        +0.j,  0.        +0.j,  0.        +0.j,
0.00249792+0.j,  0.        +0.j,  0.        +0.j,
-0.04991671+0.j,  0.        +0.j,  0.        +0.j,
-0.04991671+0.j,  0.        +0.j,  0.        +0.j,
0.99750208+0.j,  0.        +0.j,  0.        +0.j,
0.        +0.j])
"""

num_params = 1
num_wires = 4
par_domain = "R"
generator = [
qml.math.array(
[
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, -1j, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, -1j, 0, 0, -1j, 0, 0, 0, 0, 0, 0],
[0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0],
[0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1j, 0],
[0, 0, 0, 0, 0, 0, 1j, 0, 0, 1j, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1j, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1j, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
]
),
-1 / 2,
]
has_unitary_generator = False

@classmethod
def _matrix(cls, *params):
# This matrix is the "sign flipped" version of that on p18 of https://arxiv.org/abs/2104.05695,
# where the sign flip is to adjust for the opposite convention used by authors for naming wires.
# Additionally, there was a typo in the sign of a matrix element "s" at [2, 8], which is fixed here.

phi = params
c = qml.math.cos(phi / 2)
s = qml.math.sin(phi / 2)

matrix = [
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, c, 0, 0, -s, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, c, 0, 0, 0, 0, 0, -s, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, c ** 2, 0, 0, -c * s, 0, 0, -c * s, 0, 0, s ** 2, 0, 0, 0],
[0, s, 0, 0, c, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, c * s, 0, 0, c ** 2, 0, 0, -(s ** 2), 0, 0, -c * s, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, c, 0, 0, 0, 0, 0, -s, 0, 0],
[0, 0, s, 0, 0, 0, 0, 0, c, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, c * s, 0, 0, -(s ** 2), 0, 0, c ** 2, 0, 0, -c * s, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, c, 0, 0, -s, 0],
[0, 0, 0, s ** 2, 0, 0, c * s, 0, 0, c * s, 0, 0, c ** 2, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, s, 0, 0, 0, 0, 0, c, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, s, 0, 0, c, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
]

# first stack each row and then stack all the rows
U = qml.math.stack([qml.math.stack(row) for row in matrix], axis=0)

return U

[docs]    @staticmethod
def decomposition(phi, wires):
# This decomposition is the "upside down" version of that on p18 of https://arxiv.org/abs/2104.05695
decomp_ops = [
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),
qml.RY(phi / 2, wires=wires),
qml.RY(phi / 2, wires=wires),
qml.RY(phi / 2, wires=wires),
qml.RY(phi / 2, wires=wires),
qml.CNOT(wires=[wires, wires]),
qml.CNOT(wires=[wires, wires]),