Skip to main content

Documentation Index

Fetch the complete documentation index at: https://prod-mint.classiq.io/llms.txt

Use this file to discover all available pages before exploring further.

View on GitHub

Open this notebook in GitHub to run it yourself
!pip install -qq "classiq[chemistry]"


import time
from functools import reduce
from operator import mul

import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import expm

from classiq import *

​
Continuous-Time Quantum Walk on a graph

Let G(V,E)G(V, E) be a graph where VV is the vertex set and EE is the edge set, with NN representing the number of vertices (∣V∣=N|V| = N). In this model, a walker diffuses across the vertices of the graph via its edges. The distribution of the walker is represented as a wavefunction ∣ψ(t)⟩|\psi(t)\rangle in an NN-dimensional Hilbert space spanned by the basis {∣0⟩,∣1⟩,…,∣N−1⟩}\{|0\rangle, |1\rangle, \dots, |N-1\rangle\}. The probability of finding the walker at vertex vj∈Vv_j \in V is given by: pj=∣⟨j∣ψ(t)⟩∣2p_j = \left|\langle j|\psi(t)\rangle\right|^2 Given an initial state ∣ψ(0)⟩|\psi(0)\rangle, the time evolution of the continuous-time quantum walk is determined by: ∣ψ(t)⟩=e−iHt∣ψ(0)⟩|\psi(t)\rangle = e^{-iHt} |\psi(0)\rangle Here, the Hamiltonian HH is defined using the adjacency matrix AA of the graph as H=−γAH = -\gamma A, where γ>0\gamma > 0 is a constant representing the transmission rate of the walker [1].

​
Frenkel exciton Hamiltonian

\def\H{\hat{H}} \def\x{\hat{x}} \def\p{\hat{p}} \def\bin{\operatorname{bin}} \def\a{\hat{a}} \def\aD{\hat{a}^\dagger} Photosynthesis is the process by which plants use light energy to synthesize organic compounds from carbon dioxide and water. During this process, pigment-protein complexes known as light-harvesting complexes absorb light energy or receive it from other complexes, transporting this excitation energy to a reaction center. The most critical aspect of this process is the energy transfer mechanism. It is experimentally known that quantum interference plays a role in this mechanism, and one of the light-harvesting complexes where such phenomena have been observed is the Fenna-Matthews-Olson (FMO) complex [2]. The FMO complex has a trimeric structure, with each subunit containing eight Bacteriochlorophyll (BChl) a pigments. Seven of these are strongly bound within the protein scaffold and interact with each other due to their arrangement within the scaffold [3]. It is well established that the optical spectra of the FMO complex are primarily determined by the interactions within a single subunit [4]. Let us describe the movement of these excitations through Hamiltonian time evolution. Consider the ground state of the system as ∣ψ0⟩|\psi_0\rangle and the excited state localized on a single pigment (a Frenkel exciton) as ∣m⟩=\aDm∣ψ0⟩|m\rangle = \aD_m |\psi_0\rangle. In this framework, the Frenkel exciton Hamiltonian for the FMO complex is defined as follows: H=∑m=1Nεm\aDm\am+∑n<mNJmn(\aDm\an+\aDn\am),H = \sum_{m=1}^N \varepsilon_m \aD_m \a_m + \sum_{n < m}^N J_{mn} (\aD_m \a_n + \aD_n \a_m), where \aDm\aD_m and \am\a_m are the creation and annihilation operators for electronic excitation at the mm-th pigment, NN is the number of pigments, εm\varepsilon_m is the site energy of pigment mm, JmnJ_{mn} is the coupling constant between pigments mm and nn. In this notebook, we will perform a continuous-time quantum walk using the Frenkel exciton Hamiltonian given in [4]. 
fmo_hamiltonian = np.array(
    [
        [280, -106, 8, -5, 6, -8, -4],
        [-106, 420, 28, 6, 2, 13, 1],
        [8, 28, 0, -62, -1, -9, 17],
        [-5, 6, -62, 175, -70, -19, -57],
        [6, 2, -1, -70, 320, 40, -2],
        [-8, 13, -9, -19, 40, 360, 32],
        [-4, 1, 17, -57, -2, 32, 260],
    ],
    dtype=float,
)  # unit: cm^-1
The aforementioned Hamiltonian can also be written as follows: H=∑m=1Nεm∣m⟩⟨m∣+∑n<mNJmn(∣m⟩⟨n∣+∣n⟩⟨m∣)H = \sum_{m=1}^N \varepsilon_m |m\rangle\langle m| + \sum_{n < m}^N J_{mn} (|m\rangle\langle n| + |n\rangle\langle m|) To handle this Hamiltonian on a quantum computer, we must express it as a sum of Pauli strings. First, to map the qubit states to each site’s basis, we represent the integers from 00 to N−1N-1 in binary notation as follows: ∣m⟩=∣mL−1⟩⊗⋯⊗∣m1⟩⊗∣m0⟩|m\rangle = |m_{L-1}\rangle \otimes \cdots \otimes |m_1\rangle \otimes |m_0\rangle Where mi∈{0,1}m_i \in \{0, 1\}, and L=⌈log⁡2N⌉L = \lceil \log_2 N\rceil is the number of qubits required to represent the integers from 00 to N−1N-1 in binary. In this representation, the projection operators appearing in the Hamiltonian are given by: ∣m⟩⟨n∣=∣mL−1⟩⟨nL−1∣⊗⋯⊗∣m1⟩⟨n1∣⊗∣m0⟩⟨n0∣|m\rangle\langle n| = |m_{L-1}\rangle\langle n_{L-1}| \otimes \cdots \otimes |m_1\rangle\langle n_1| \otimes |m_0\rangle\langle n_0| Consequently, the Hamiltonian can be expressed as a sum of Pauli strings using the following transformations: ∣0⟩⟨0∣=12(I+Z),∣1⟩⟨0∣=12(X−iY),∣1⟩⟨1∣=12(I−Z),∣0⟩⟨1∣=12(X+iY).\begin{equation*} \begin{split} |0\rangle\langle 0| & = \frac{1}{2}(I + Z),\\ |1\rangle\langle 0| & = \frac{1}{2}(X - iY),\\ \end{split} \qquad \begin{split} |1\rangle\langle 1| & = \frac{1}{2}(I - Z),\\ |0\rangle\langle 1| & = \frac{1}{2}(X + iY).\\ \end{split} \end{equation*} Note that the Frenkel exciton Hamiltonian restricts the excitation to the single-exciton manifold, so it is not scalable. However, for a better description of energy-transfer dynamics, we need to mix single-exciton states with ground state or multiexciton states, and in this case, the computational complexity scales exponentially with the number of pigments [5]. 
def create_projector(ket: int, bra: int, qubits: list[int]) -> SparsePauliOp:
    """Create projector operator for given ket and bra states.

    Args:
        ket (int): Ket state.
        bra (int): Bra state.
        qubits (list[int]): List of qubit indices.

    Returns:
        SparsePauliOp: Projector operator.
    """
    m = len(qubits)
    if not (0 <= ket < 2**m and 0 <= bra < 2**m):
        raise ValueError("Ket and bra states must be in the range [0, 2^m).")
    ket_bin = bin(ket)[2:].zfill(m)
    bra_bin = bin(bra)[2:].zfill(m)
    paulis_coeffs = []
    for i in range(m):
        # little endian convention:
        # the least significant bit corresponds to the first qubit in the list
        if ket_bin[m - i - 1] == "0" and bra_bin[m - i - 1] == "0":
            paulis_coeffs.append((Pauli.I(qubits[i]) + Pauli.Z(qubits[i])) / 2)
        elif ket_bin[m - i - 1] == "0" and bra_bin[m - i - 1] == "1":
            paulis_coeffs.append((Pauli.X(qubits[i]) + 1j * Pauli.Y(qubits[i])) / 2)
        elif ket_bin[m - i - 1] == "1" and bra_bin[m - i - 1] == "0":
            paulis_coeffs.append((Pauli.X(qubits[i]) - 1j * Pauli.Y(qubits[i])) / 2)
        elif ket_bin[m - i - 1] == "1" and bra_bin[m - i - 1] == "1":
            paulis_coeffs.append((Pauli.I(qubits[i]) 

- Pauli.Z(qubits[i])) / 2)

    return reduce(mul, paulis_coeffs)


fmo_hamiltonian_sparse = []
for i in range(fmo_hamiltonian.shape[0]):
    for j in range(fmo_hamiltonian.shape[1]):
        fmo_hamiltonian_sparse.append(
            create_projector(i, j, list(range(3))) * fmo_hamiltonian[i, j]
        )
fmo_hamiltonian_sparse = reduce(lambda x, y: x + y, fmo_hamiltonian_sparse)
To approximate the time evolution of the constructed Hamiltonian, we employ the first-order Suzuki-Trotter decomposition: e−itH=exp⁡{−it∑j=1NhjHj}≈(∏jNe−ithjHj/r)r.\begin{equation*} e^{-itH}=\exp\left\{-it\sum_{j=1}^N h_j H_j\right\} \approx \left(\prod_{j}^N e^{-it h_j H_j/r}\right)^r. \end{equation*} Here, we fix the number of time slices r=10r = 10. 
@qfunc
def main(qnum: Output[QNum], time_point: CReal):
    allocate(3, qnum)
    suzuki_trotter(
        fmo_hamiltonian_sparse,
        evolution_coefficient=time_point,
        order=1,
        repetitions=10,
        qbv=qnum,
    )


qprog = synthesize(main)
qprog = set_quantum_program_execution_preferences(
    qprog,
    preferences=ExecutionPreferences(
        num_shots=1,
        backend_preferences=ClassiqBackendPreferences(
            backend_name="simulator_statevector"
        ),
    ),
)
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/3B5kTnagDW5FSbG1FiSEfR58Zlw
  


time_points = np.linspace(0, 0.1, 20)
trotter_results = []
durations = []

with ExecutionSession(qprog) as es:
    for time_point in time_points:
        time1 = time.time()

        result = es.sample({"time_point": time_point})

        probs = [0.0] * 2**3

        for detail in result.parsed_state_vector:
            probs[detail.state["qnum"]] = abs(detail.amplitude) ** 2

        trotter_results.append(probs)
        time2 = time.time()
        duration = time2 - time1
        durations.append(duration)
        print(f"Time point: {time_point:.4f}, Duration: {duration:.4f} seconds")
Output:

Time point: 0.0000, Duration: 1.7441 seconds
  Time point: 0.0053, Duration: 1.3106 seconds
  Time point: 0.0105, Duration: 1.3255 seconds
  Time point: 0.0158, Duration: 1.2085 seconds
  Time point: 0.0211, Duration: 0.9949 seconds
  Time point: 0.0263, Duration: 0.8305 seconds
  Time point: 0.0316, Duration: 1.4334 seconds
  Time point: 0.0368, Duration: 1.3218 seconds
  Time point: 0.0421, Duration: 0.9502 seconds
  Time point: 0.0474, Duration: 1.1008 seconds
  Time point: 0.0526, Duration: 1.5134 seconds
  Time point: 0.0579, Duration: 1.3479 seconds
  Time point: 0.0632, Duration: 1.3225 seconds
  Time point: 0.0684, Duration: 1.2989 seconds
  Time point: 0.0737, Duration: 0.8707 seconds
  Time point: 0.0789, Duration: 1.2082 seconds
  Time point: 0.0842, Duration: 0.9853 seconds
  Time point: 0.0895, Duration: 1.0298 seconds
  Time point: 0.0947, Duration: 1.4580 seconds
  Time point: 0.1000, Duration: 1.3582 seconds
Let’s compare the calculation results with those obtained using scipy. Here, we plot the time evolution of the exciton population for pigments 0, 1, and
We see that accuracy decreases over time due to decomposition error. 
def simulate_ctqw(time_steps, initial_node=0):
    H = fmo_hamiltonian

    # initial state
    psi0 = np.zeros(fmo_hamiltonian.shape[0], dtype=complex)
    psi0[initial_node] = 1.0

    results = []
    for t in time_steps:
        # time evolution operator: U = exp(-i * H * t)
        U = expm(-1j * H * t)
        psi_t = np.dot(U, psi0)

        probability = np.abs(psi_t) ** 2
        results.append(probability.tolist())

    return results


t_max = 0.1
time_steps = np.linspace(0, t_max, 100)

probs = simulate_ctqw(time_steps, initial_node=0)

plt.plot(
    time_points, [p[0] for p in trotter_results], "bs--", label="Pigment 0 (Classiq)"
)
plt.plot(
    time_points, [p[1] for p in trotter_results], "g^--", label="Pigment 1 (Classiq)"
)
plt.plot(
    time_points, [p[2] for p in trotter_results], "rx--", label="Pigment 2 (Classiq)"
)

plt.plot(time_steps, [p[0] for p in probs], "b-", label="Pigment 0 (Scipy)")
plt.plot(time_steps, [p[1] for p in probs], "g-", label="Pigment 1 (Scipy)")
plt.plot(time_steps, [p[2] for p in probs], "r-", label="Pigment 2 (Scipy)")
plt.title("Continuous-Time Quantum Walk with FMO Hamiltonian")
plt.xlabel("Time")
plt.ylabel("Population")
plt.legend(bbox_to_anchor=(1.02, 1))
plt.grid(True, linestyle="--", alpha=0.7)
plt.show()
output

​
References

[1]: R. Portugal and J. K. Moqadam, Efficient circuit implementations of continuous-time quantum walks for quantum search, Entropy (Basel) 27, 454 (2025). [2]: M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, Environment-assisted quantum walks in photosynthetic energy transfer, J. Chem. Phys. 129, 174106 (2008). [3]: M. Maiuri, E. E. Ostroumov, R. G. Saer, R. E. Blankenship, and G. D. Scholes, Coherent wavepackets in the Fenna-Matthews-Olson complex are robust to excitonic-structure perturbations caused by mutagenesis, Nat. Chem. 10, 177 (2018). [4]: M. Cho, H. M. Vaswani, T. Brixner, J. Stenger, and G. R. Fleming, Exciton analysis in 2D electronic spectroscopy, J. Phys. Chem. B 109, 10542 (2005). [5]: C.-K. Lee, J. W. Zhong Lau, L. Shi, and L. C. Kwek, Simulating energy transfer in molecular systems with digital quantum computers, J. Chem. Theory Comput. 18, 1347 (2022).