June 6th, 2025

Introduction to PyOR: Simulating a Spectrum

Author: Vineeth Francis Thalakottoor Jose Chacko
Editor: Quentin Stern

Introduction

Welcome to the second post on PyOR!

In the previous post, I introduced some of PyOR’s key features and its approach to representing quantum objects. In this post, I want to give an overview of PyOR and walk through a concrete example: how to simulate a spectrum using PyOR.

We’ll consider a simple two-spin system with J-coupling and go through the simulation process step by step. This includes:

Import Packages (PyOR and Python)
Creating a spin system
Generating spin operators
Building the Hamiltonian and density matrix
Applying a pulse
Evolving the system
Plotting the final spectrum

First, I recommend looking at the complete script, provided in the form of a Jupyter Notebook.

To download this script, click here. To run the simulation, you’ll need the PyOR package, which includes two folders: Source_Doc and Examples. In the script, make sure to set the path to the Source_Doc folder on your local machine.

For a smooth experience, I recommend using Anaconda and Visual Studio Code. There are many tutorials available online on how to connect Anaconda and Jupyter Notebook to Visual Studio Code. Also, make sure to install ipympl (pip install ipympl) for interactive Matplotlib features.

Import PyOR modules and other Python packages

First, set the variable SourcePath by specifying the path to the Source_Doc directory that contains all the modules of the PyOR package.

# Define the source path
SourcePath = 'Path to /Source_Doc'

# Import Python packages and add source path
import sys
sys.path.append(SourcePath)
import time
import numpy as np
%matplotlib ipympl

# Import PyOR package
from PyOR_QuantumSystem import QuantumSystem as QunS
from PyOR_HardPulse import HardPulse
from PyOR_Evolution import Evolutions
from PyOR_Plotting import Plotting
import PyOR_SignalProcessing as Spro

import sys: Imports Python’s sys module, which provides access to system-specific parameters and functions.
sys.path: This is a list of directories where Python looks for modules to import. By default, it includes the current directory, standard library paths, and any installed packages.
sys.path.append(SourcePath): Adds your custom directory (stored in the variable SourcePath) to the list of module search paths.
import time: Imports Python’s built-in time module.
import numpy as np: Imports the NumPy library and assigns it the alias np.
%matplotlib ipympl: A Jupyter Notebook magic command. Enables interactive plotting using ipympl as the Matplotlib backend.

In this simulation, we use the following PyOR modules (Python files containing either classes or functions): PyOR_QuantumSystem (contain simulation parameters set by user, create spin operators and quantum states), PyOR_HardPulse (to apply pulse), PyOR_Evolution (for evolving the system), PyOR_SignalProcessing (for the Fourier transform), and PyOR_Plotting (for data visualization).

Define the spin system

At the core of PyOR lies the powerful QuantumSystem class, defined in the PyOR_QuantumSystem module. This class is responsible for generating all the necessary spin operators used in magnetic resonance simulations. The spin operators include three components of angular momentum operators (say $I_x$ , $I_y$ and $I_z$ ) and two ladder operators ( $I_p = I_x + j I_y$ and $I_m = I_x - j I_y$ where $j = \sqrt{-1}$ .

The first step is to define your spin system. For example, if you’re working with two protons, you might label them “A” and “B”. This is done using a simple Python dictionary:

Spin_list = {“A” : “H1”, “B” : “H1”}

This Spin_list, a python dictionary is the input to the QuantumSystem class, typically imported under the alias QunS. You can then instantiate your quantum system like this:

# Define the spin system
Spin_list = {"A" : "H1", "B" : "H1"}

QS = QunS(Spin_list,PrintDefault=False)

# initialize the system
QS.Initialize()

Once this is done, QS contains all the spin operators and internal mappings needed to begin your magnetic resonance simulations. The method Initialize() calls all other methods responsible for creating the spin operators:

SpinOperatorsSingleSpin: Create the angular momentum spin operator for a single spin with an arbitrary spin quantum number, which will be used by the methods below.
SpinOperator: Create all six spin operators mentioned previously.
SpinOperator_Sub: Create all six spin operators for the subsystem.

Spin Operators

The initialized spin operators can be accessed as follows. For example, the $z$ component of the spin operator for the spin labeled “A” can be called using:

# Call spin operator of the system, example: z spin operator of spin 'A'
QS.Az
# Output: <PyOR_QuantumObject.QunObj at 0x7f5b6494e1e0>

QS.Az.matrix

# Call spin operator of the sub system, example: z spin operator of spin 'A'
QS.Az_sub
# Output: <PyOR_QuantumObject.QunObj at 0x7f5b5b801d30>

QS.Az_sub.matrix

QS.Az is a quantum object as mentioned in my previous post, which represents the spin operator $A_z$ . For symbolic representation, QS.Az.matrix.

$\left[ \begin{array}{cccc} 0.5 & 0 & 0 & 0 \\ 0 & 0.5 & 0 & 0 \\ 0 & 0 & -0.5 & 0 \\ 0 & 0 & 0 & -0.5 \end{array} \right]$

PyOR not only creates spin operators for the entire system, but also for individual subsystems. To view the symbolic representation of the $z$ component of the spin operator in the subspace of spin “A”, use QS.Az_sub.matrix.

$\left[ \begin{array}{cc} 0.5 & 0 \\ 0 & -0.5 \end{array} \right]$

Set the simulation parameters

After creating the spin operators, the next step is to set all the simulation parameters (as the attributes of the QuantumSystem class), such as the vector space, master equation, $B_0$ field, chemical shift values, J-coupling constants, relaxation parameters, among others. We will explore these in upcoming posts. For now, we will begin with the basic parameters. We are going to simulate two spin system (protons) that are J-coupled.

# Master Equation
QS.PropagationSpace = "Hilbert"
QS.MasterEquation = "Redfield"

# Operator Basis
QS.Basis_SpinOperators_Hilbert = "Zeeman"

# Relaxation Process
QS.Rprocess = "Phenomenological"
QS.R1 = 1
QS.R2 = 1

QS.Update()

In this post, I introduce four parameters:

QS.PropagationSpace: Choose the vector space, either “Schrodinger” or “Hilbert” or “Liouville ” space.
QS.MasterEquation: Choose the master equation, either “Redfield” or “Lindblad” equation.
QS.Basis_SpinOperators_Hilbert: Choose the basis of the spin operators, either “Zeeman” or “Singlet Triplet”.
QS.Rprocess: The relaxation process. We will see “Phenomenological” relaxation this time. In future posts, we will see other possible relaxations implemented in PyOR. Here QS.R1 is the longitudinal relaxation rate and QS.R2 the relaxation rate for all the coherence, which is set to 1 Hz.

Generate the system Hamiltonian

In PyOR, we can create different types of Hamiltonians like Zeeman, J-coupling, Dipolar, chemical shift anisotropy (CSA), Quadrupole, among others, which are defined in the PyOR_Hamiltonian module. But will see how to create Hamiltonians (Zeeman and J-coupling) from spin operators, see below:

# Zeeman Hamiltonian (rotating frame)
Hz = 2.0 * np.pi * 10 * QS.Az + 2.0 * np.pi * 50 * QS.Bz

# J coupling Hamiltonian
Hj = 2.0 * np.pi * 5 * (QS.Ax * QS.Bx + QS.Ay * QS.By + QS.Az * QS.Bz)

Zeeman Hamiltonian

The interaction of spins with the external static field is given by:

$H_z = \omega_A \cdot I_{z}^{A} + \omega_B \cdot I_{z}^{B}$

where $\omega_A = 2\pi \cdot f_A$ , $f_A$ is the Larmor frequency of spin A, similarly for spin B. In the script above,

Hz = 2.0 * np.pi * 10 * QS.Az + 2.0 * np.pi * 50 * QS.Bz

generate the Zeeman Hamiltonian of spins A and B with Larmor frequencies of 10 Hz and 50 Hz, respectively.

J coupling Hamiltonian

The through-bond scalar (J) coupling between spins A and B is given by:

$H_J = 2\pi \cdot J \cdot \left( I_{x}^{A} I_{x}^{B} + I_{y}^{A} I_{y}^{B} + I_{z}^{A} I_{z}^{B} \right)$

where $J$ is the J-coupling constant. To generate a J coupling Hamiltonian with a coupling constant of 5 Hz:

Hj = 2.0 * np.pi * 5 * (QS.Ax * QS.Bx + QS.Ay * QS.By + QS.Az * QS.Bz)

Generate the equilibrium density matrix

Liouville–von Neumann equation

In Nuclear Magnetic Resonance (NMR), the state of the system is described by a density matrix. The diagonal elements of the density matrix represent the populations of each quantum state, while the off-diagonal elements capture the coherence or phase relationships between different quantum states. The evolution of the density matrix provides insight into how the system changes over time under the influence of a given Hamiltonian. The following equation (Liouville–von Neumann equation or in short, Liouville equation) governs this evolution:

$\frac{d}{dt}\rho = -j [H, \rho]$

However, the Liouville–von Neumann equation alone does not account for the system’s relaxation toward thermal equilibrium. To describe relaxation processes, an additional term must be included in the equation to model the interaction between the system and its environment. A commonly used approach in NMR to incorporate relaxation is the Redfield master equation, which is shown below:

$\frac{d}{dt} \rho = -j [H, \rho] + R(\rho - \rho_{eq})$

Here, the second part $R$ corresponds to the relaxation matrix and is a function of $\rho - \rho_{eq}$ . $\rho_{\text{eq}}$ is the equilibrium density matrix. In this work, we use a simple relaxation model based on a phenomenological relaxation matrix, $R$ . The diagonal elements of this matrix are given by QS.R1, representing longitudinal relaxation, while the off-diagonal elements are QS.R2, corresponding to transverse relaxation (also other coherences). Thus, the relaxation part in the Redfield master equation can be simplified to the relaxation matrix R time $\rho - \rho_{eq}$

$\frac{d}{dt}\rho = -j [H, \rho] + R * (\rho - \rho_{eq})$

Equilibrium density matrix

We consider that our system is in an equilibrium state, defined by QS.Az + QS.Bz, let this be our initial state, rho_in. After that, we apply a hard pulse on the spins and make them out of equilibrium. In the simulation, we define an initial (rho_in) and a final state (rhoeq) where the system remains in equilibrium. See below code snippet:

# Initial Density Matrix
rho_in = QS.Az + QS.Bz
rho_in.matrix

# Final Density Matrix
rhoeq = QS.Az + QS.Bz
rhoeq.matrix

Applying the hard pulse

After defining an initial state, let’s apply a 90-degree pulse. To do this, we need to import the HardPulse class and create an instance HardP. The method, Rotate_Pulse takes arguments: density matrix, the flip angle, and the operator about which the user needs to make the rotation. In this simulation, we choose to rotate the spins about the y-axis, ie. about the operator QS.Ay + QS.By.

Note: Numerical computation comes with risk of numerical errors, so to mitigate it, the method Tolarence of QunObj (which I introduced in the last post) can be useful. Tolarence(1.0e-5) can make any small matrix element into zero by defining a tolerance value as an argument, here it is 1.0e-5.

HardP = HardPulse(QS)

flip_angle = 90.0   # Flip angle
rho = HardP.Rotate_Pulse(rho_in,flip_angle,QS.Ay + QS.By).Tolarence(1.0e-5)
rho.matrix

Evolution of the density matrix

Here comes the most important part of this post, evolution. For this, you need to make an instance of the Evolution class, EVol. In PyOR you can solve the Liouville equation either by using a unitary propagator or by solving the set of ordinary differential equations (OEDs). In Hilbert space, the Redfield or Lindblad master equation is supported only by ODE solver and not by a unitary propagator. PyOR uses Scipy solve_ivp for solving a set of differential equations.

The user can define the dwell time, Acquisition time, solver method, and propagation by setting the QS.AcqDT, QS.AcqAQ, QS.OdeMethod and QS.PropagationMethod parameters respectively.

QS.PropagationMethod comes with a lot of options, a few are mentioned here, and we will see them in later posts:

“Unitary Propagator”
“Unitary Propagator Time Dependent”
“ODE Solver”
“ODE Solver Lindblad”

The method Evolution propagate the density matrix with arguments density matrix to propagate, the equilibrium density matrix, and the Hamiltonian, defined in the Evolution class. It returns the array of time points and the density matrix for each time point.

t, rho_t = EVol.Evolution(rho,rhoeq,Hz+Hj)

QS.AcqDT = 0.0001
QS.AcqAQ = 5.0
QS.OdeMethod = 'DOP853'
QS.PropagationMethod = "ODE Solver"

EVol = Evolutions(QS)

start_time = time.time()
t, rho_t = EVol.Evolution(rho,rhoeq,Hz+Hj)
end_time = time.time()
timetaken = end_time - start_time
print("Total time = %s seconds " % (timetaken))

Expectation values

The expectation value of any observable, say, QS.Ap + QS.Bp (to observe Mt, transverse magnetization) or QS.Az + QS.Bz (to observe Mz, longitudinal magnetization) can be found using the method Expectation (argument: array of density matrix and observable), which returns the expectation value along with the time points.

det_Mt = QS.Ap + QS.Bp
det_Z = QS.Az + QS.Bz

t, Mt = EVol.Expectation(rho_t,det_Mt)
t, Mz = EVol.Expectation(rho_t,det_Z)

Plotting the transverse and longitudinal magnetization

PyOR has many methods to make plots for easy visualization, defined in the Plotting class. Let’s see one method:

Plotting_SpanSelector(t,Mt,"time (s)","Mt","red")

where the arguments: t for the x-axis data, Mt for the y-axis data, "time (s)" as the x-axis label, "Mt" as the y-axis label, and "red" as the color for the plot.

plot = Plotting(QS)

plot.Plotting_SpanSelector(t,Mt,"time (s)","Mt","red",saveplt=True,savename="Transverse Magnetization")

plot.Plotting_SpanSelector(t,Mz,"time (s)","Mz","red",saveplt=True,savename="Longitudinal Magnetization")

Transverse magnetization

Longitudinal magnetization

Plotting the spectra

The FourierTransform function in the PyOR_SignalProcessing module can be used to calculate the Fourier transform of the transverse magnetization.

freq, spectrum = Spro.FourierTransform(Mt,QS.AcqFS,N)

where the arguments are the data to be Fourier transformed, QS.AcqFS = 1/QS.AcqDT, and zero filling (that means the length of the transformed axis is N times the length of the data)

freq, spectrum = Spro.FourierTransform(Mt,QS.AcqFS,5)

plot.Plotting_SpanSelector(freq,spectrum,"Frequency (Hz)","Spectrum","red",saveplt=True,savename="Spectrum")

Spectra

Conclusion

In this post, we saw how to define a spin system, generate spin operators, create the Hamiltonians, apply a hard pulse, evolve the density matrix, calculate the expectation, and perform data visualization, using PyOR. In coming posts, we will see different Hamiltonians that are defined in the Hamiltonian class.

// Initialize the counter variable
$counter = 1;

// Loop from 1 to 10
while ($counter <= 10) {
// Print the current value of the counter
echo $counter . "\n";

// Increment the counter by 1
$counter++;
}
?>

1 2 3 4 5 6 7 8 9 10

#include

int main() {
// Initialize the counter variable
int counter = 1;

// Loop from 1 to 10
while (counter <= 10) {
// Print the current value of the counter
printf("%d\n", counter);

// Increment the counter by 1
counter++;
}

return 0;
}

1 2 3 4 5 6 7 8 9 10

#include

int main() {
// Initialize the counter variable
int counter = 1;

// Loop from 1 to 10
while (counter <= 10) {
// Print the current value of the counter
std::cout << counter << std::endl;

// Increment the counter by 1
counter++;
}

return 0;
}

1 2 3 4 5 6 7 8 9 10

// Initialize the counter variable
let counter = 1;

// Loop from 1 to 10
while (counter <= 10) {
// Print the current value of the counter
console.log(counter);

// Increment the counter by 1
counter++;
}

Introduction to PyOR: Simulating a Spectrum

Introduction

Import PyOR modules and other Python packages

Define the spin system

Spin Operators

Set the simulation parameters

Generate the system Hamiltonian

Zeeman Hamiltonian

J coupling Hamiltonian

Generate the equilibrium density matrix

Liouville–von Neumann equation

Further reading

Equilibrium density matrix

Applying the hard pulse

Evolution of the density matrix

Expectation values

Plotting the transverse and longitudinal magnetization

Transverse magnetization

Longitudinal magnetization

Plotting the spectra

Spectra

Conclusion

Leave a Reply Cancel reply

Artist Credit:

Contact:

Website designed and developed by:

Introduction to PyOR: Simulating a Spectrum

Introduction

Import PyOR modules and other Python packages

Define the spin system

Spin Operators

Set the simulation parameters

Generate the system Hamiltonian

Zeeman Hamiltonian

J coupling Hamiltonian

Generate the equilibrium density matrix

Liouville–von Neumann equation

Further reading

Equilibrium density matrix

Applying the hard pulse

Evolution of the density matrix

Expectation values

Plotting the transverse and longitudinal magnetization

Transverse magnetization

Longitudinal magnetization

Plotting the spectra

Spectra

Conclusion

Leave a Reply Cancel reply

Artist Credit:

Contact:

Website designed and developed by:

Quantum Insights