# Converting circuit gates to Clifford + T gates

Dear pennylane team,
I have created a circuit by performing the following operations:

First, I created the Hamiltonian.

``````symbols = ["H", "O", "H"]
coordinates = np.array([-0.0399, -0.0038, 0.0, 1.5780, 0.8540, 0.0, 2.7909, -0.5159, 0.0])
hamiltonian, num_qubits = qchem.molecular_hamiltonian(symbols, coordinates)
``````

Second, I passed the Hamiltonian down to the `ApproxTimeEvolution` function.

``````wires = range(num_qubits)
num_wires = num_qubits

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

@qml.qnode(dev)
def circuit(time):
ApproxTimeEvolution(hamiltonian, time, 1)
return [qml.expval(qml.PauliZ(wires=i)) for i in wires]
``````

Next, I would like to see what the circuit is made of. So I run the `spec_func` and get the following not-informative output:

``````{'gate_sizes': defaultdict(int, {14: 1}),
'gate_types': defaultdict(int, {'ApproxTimeEvolution': 1}),
'num_operations': 1,
'num_observables': 14,
'num_diagonalizing_gates': 0,
'num_used_wires': 14,
'depth': 1,
'num_trainable_params': 1,
'num_device_wires': 14,
'diff_method': 'best',
``````

The problem with this is that it doesn’t show what gate types are used but just says `ApproxTimeEvolution`.

So, I have two questions:
1. How do I actually see what kind of gates are used in the circuit?
2. I would like the circuit to use exclusively the Clifford + T gate set. How do I do this?

Thanks!

Hi @Einar_Gabbassov, it’s great to see you in the Forum!

1. You can add the keyword argument `expansion_strategy='device'`.
You can also use `print(qml.draw(circuit, expansion_strategy='device')(time))` and replace `time` with the value you want.

2. If you want your circuit to use exclusively a set of gates then you can find a device which uses such gates or create a custom device with the gates that you want. You can learn how to create your own device in this blog post. Then `expansion_strategy='device'` should decompose the circuit into the gates in that device.

Please let me know if this works for you!

Hi @CatalinaAlbornoz, thanks for tips.

Unfortunattely, the suggested workflow doesn’t yield the desired gates.
It actually produces something else. Here is the code:

``````import numpy as np
import pennylane as qml
from pennylane import qchem
from pennylane.templates import ApproxTimeEvolution

class MyDevice(qml.QubitDevice):
name = 'Clifford_T_device'
short_name = "custom.qubit"
pennylane_requires = ">=0.23"
version = "0.0.1"
author = "Maurice"
# Here, I define the desired gate set to be used
operations = {"H", "S", "CNOT", "T"}
observables = {"PauliX", "PauliZ"}

def __init__(self, num_wires):
super().__init__(wires=num_wires, shots=None)

# create the initial state
self._state = np.array([1, 0])

# create a variable for future copies of the state
self._pre_rotated_state = None

@property
def state(self):
return self._pre_rotated_state

@classmethod
def capabilities(cls):
capabilities = super().capabilities().copy()
capabilities.update(
returns_state=True,
supports_finite_shots=False,
supports_tensor_observables=False
)
return capabilities

def apply(self, operations, rotations=None, **kwargs):
for op in operations:
# We update the state by applying the matrix representation of the gate
self._state = qml.matrix(op) @ self._state

# store the pre-rotated state
self._pre_rotated_state = self._state.copy()

# apply the circuit rotations
for rot in rotations or []:
self._state = qml.matrix(rot) @ self._state

def analytic_probability(self, wires=None):
if self._state is None:
return None

real = self._real(self._state)
imag = self._imag(self._state)
prob = self.marginal_prob(real ** 2 + imag ** 2, wires)
return prob

def reset(self):
"""Reset the device"""
self._state = np.array([1, 0])

hamiltonian, num_qubits = qchem.molecular_hamiltonian(symbols, coordinates)
wires = range(num_qubits)
num_wires = num_qubits

dev = MyDevice(num_wires)

time=1
trotter_steps=1
@qml.qnode(dev)
def circuit(time):
ApproxTimeEvolution(hamiltonian, time, trotter_steps)
return [qml.expval(qml.PauliZ(wires=i)) for i in wires]

specs_func = qml.specs(circuit, expansion_strategy='device')
print(specs_func(1))
t = circuit.qtape.to_openqasm()
print(t)

``````

Even though the indicated gate set is universal, I feel like for the suggested approach to work we need some sort of Solovay-Kitaev algorithm, which will assist the conversion process. I’m not sure if the suggested approach anyhow takes this into account. Please correct me if I’m wrong.

For you convenience, I also attach the output

``````OPENQASM 2.0;
include "qelib1.inc";
qreg q[4];
creg c[4];
rz(0.4740443816356016) q[0];
rz(0.4740443816356016) q[1];
cx q[1],q[0];
rz(0.36906324698635556) q[0];
cx q[1],q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rx(-1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rx(-1.5707963267948966) q[3];
rx(1.5707963267948966) q[0];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(-0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rx(-1.5707963267948966) q[0];
rx(-1.5707963267948966) q[1];
rz(1.5707963267948966) q[2];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[2];
rx(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(-0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rz(1.5707963267948966) q[1];
rx(1.5707963267948966) q[1];
rz(1.5707963267948966) q[1];
rx(-1.5707963267948966) q[2];
rx(-1.5707963267948966) q[3];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[1];
rx(1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
cx q[3],q[2];
cx q[2],q[1];
cx q[1],q[0];
rz(0.08210294405757057) q[0];
cx q[1],q[0];
cx q[2],q[1];
cx q[3],q[2];
rz(1.5707963267948966) q[0];
rx(1.5707963267948966) q[0];
rz(1.5707963267948966) q[0];
rx(-1.5707963267948966) q[1];
rx(-1.5707963267948966) q[2];
rz(1.5707963267948966) q[3];
rx(1.5707963267948966) q[3];
rz(1.5707963267948966) q[3];
rz(-0.9211542548980106) q[2];
cx q[2],q[0];
rz(0.2812358847255707) q[0];
cx q[2],q[0];
rz(-0.9211542548980107) q[3];
cx q[3],q[0];
rz(0.3633388287831413) q[0];
cx q[3],q[0];
cx q[2],q[1];
rz(0.3633388287831413) q[1];
cx q[2],q[1];
cx q[3],q[1];
rz(0.2812358847255707) q[1];
cx q[3],q[1];
cx q[3],q[2];
rz(0.38350121257667436) q[2];
cx q[3],q[2];
measure q[0] -> c[0];
measure q[1] -> c[1];
measure q[2] -> c[2];
measure q[3] -> c[3];
``````

Hey @Einar_Gabbassov,

Does the `decomposition()` method work for you? E.g.,

``````>>> hamiltonian, num_qubits = qchem.molecular_hamiltonian(["H", "H"], np.array([[0, 0, 0.66], [0, 0, -0.66]]))
>>> qml.ApproxTimeEvolution(hamiltonian, 1, 1).decomposition()
``````

If those operations aren’t decomposed enough, you can decompose them further:

``````>>> for op in qml.ApproxTimeEvolution(hamiltonian, 1, 1).decomposition():
...     print(op.decomposition())
...
The call `qml.ApproxTimeEvolution(hamiltonian, 1, 1).decomposition()` outputs the list of gates. This is great! Thanks. Also, thanks a lot for the tip on how to decompoe even further. This was extremly useful.