Brief Example

We now provide a brief example to demonstrate the similarity between QuantumToolbox.jl and QuTiP.

CPU Computation

Let's consider a quantum harmonic oscillator with a Hamiltonian given by:

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}$$

where $\hat{a}$ and ${\hat{a}}^{†}$ are the annihilation and creation operators, respectively. We can define the Hamiltonian as follows:

using QuantumToolbox

N = 20 # cutoff of the Hilbert space dimension
ω = 1.0 # frequency of the harmonic oscillator

a = destroy(N) # annihilation operator

H = ω * a' * a

We now introduce some losses in a thermal environment, described by the Lindblad master equation:

$$\frac{d\hat{\rho }}{dt}=-i[\hat{H},\hat{\rho }]+\gamma \mathcal{D}[\hat{a}]\hat{\rho }$$

where $\hat{\rho }$ is the density matrix, $\gamma$ is the damping rate, and $\mathcal{D}[\hat{a}]$ is the Lindblad dissipator, defined as:

$$\mathcal{D}[\hat{a}]\hat{\rho }=\hat{a}\hat{\rho }{\hat{a}}^{†}-\frac{1}{2}{\hat{a}}^{†}\hat{a}\hat{\rho }-\frac{1}{2}\hat{\rho }{\hat{a}}^{†}\hat{a}$$

Lindblad master equation

See here for more details about Lindblad master equation.

We now compute the time evolution of the system using the mesolve function, starting from the initial state $|\psi (0)⟩=|3⟩$:

γ = 0.1 # damping rate

ψ0 = fock(N, 3) # initial state

tlist = range(0, 10, 100) # time list

c_ops = [sqrt(γ) * a]
e_ops = [a' * a]

sol = mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops)

We can extract the expectation value of the number operator ${\hat{a}}^{†}\hat{a}$ with the command sol.expect, and the states with the command sol.states.

GPU Computation

Extension for CUDA.jl

QuantumToolbox.jl provides an extension to support GPU computation. To trigger the extension, you need to install and import CUDA.jl together with QuantumToolbox.jl. See here for more details.

using QuantumToolbox
using CUDA
CUDA.allowscalar(false) # Avoid unexpected scalar indexing

We can easily pass the computation to the GPU, by simply passing all the QuantumObjects to the GPU:

a_gpu = cu(destroy(N)) # The only difference in the code is the cu() function

H_gpu = ω * a_gpu' * a_gpu

ψ0_gpu = cu(fock(N, 3))

c_ops = [sqrt(γ) * a_gpu]
e_ops = [a_gpu' * a_gpu]

sol = mesolve(H_gpu, ψ0_gpu, tlist, c_ops, e_ops = e_ops)