Calculating NMR Spectra with SimpCalc¶
This tutorial demonstrates how to use SimpCalc to design and run NMR simulations using SIMPSON. We will cover defining spin systems, setting up the calculator, using different pulse sequences, and handling simulation outputs.
1. Defining the Spin System¶
The spin system defines the nuclei, their chemical shifts, and couplings. You can define this using a SIMPSON-formatted string or use soprano objects if you have them from DFT calculations.
from __future__ import annotations
import shutil
import matplotlib.pyplot as plt
from simpyson.calculator import SimpCalc
# Check if SIMPSON is installed
SIMPSON_INSTALLED = shutil.which('simpson') is not None
if not SIMPSON_INSTALLED:
print("WARNING: 'simpson' executable not found in PATH. You can generate input files but cannot run simulations directly.")
else:
print("SIMPSON detected. Ready to run simulations.")
2. Setting up the Calculator¶
SimpCalc is the main interface. You initialize it with your spin system and simulation parameters.
Note: SimpCalc requires several parameters to be explicitly set, even if standard defaults are desired.
sim = SimpCalc(
spinsys=spinsys_str,
proton_frequency='400MHz',
sw=2000, # Spectral width in Hz
np=4096, # Number of points
variable_ref=0, # Center frequency
# Creating standard static simulation parameters
spin_rate='10e3', # MAS rate in Hz
start_operator="I1x",
detect_operator="I1p",
method="direct",
crystal_file="rep100", # Powder averaging
gamma_angles=10,
verbose=0
)
# Save file
sim.save('../../examples/calculator/my_simulation.in')
3. Pulse Sequences¶
SimpCalc supports various pulse sequences. The default is 'no_pulse' (just acquire). You can specify others like 'pulse90' or 'cpmas'.
3.1 Standard Pulse-Acquire (Pulse90)¶
This uses a simple 90-degree pulse before acquisition.
from simpyson.templates import Pulse90
pulse90 = Pulse90()
spinsys_13C = """
channels 13C
nuclei 13C 13C
shift 1 50p 0 0 0 0 0
shift 2 20p 0 0 0 0 0
"""
sim_90 = SimpCalc(
spinsys=spinsys_13C,
proton_frequency='400e6',
sw=20000,
np=4096,
# Standard Parameters
spin_rate='10e3', # MAS rate in Hz
start_operator="Inz", # Start with z, then pulse90 flips it
detect_operator="Inp",
method="direct",
crystal_file="rep100",
gamma_angles=10,
verbose=0,
pulse_sequence=pulse90,
lb=10, # Add line broadening
zerofill=4096 # Add zerofill
)
4. Running the Simulation¶
If simpson is available, you can run the simulation explicitly. The run method returns a Simpy object containing the result.
output = sim_90.run()
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(15,5), constrained_layout=True)
ax[0].plot(output.fid['time'], output.fid['real'])
ax[0].set_title('FID')
ax[0].set_xlabel('Time (ms)')
ax[1].plot(output.spe['ppm'], output.spe['real'])
ax[1].set_title('Spectrum in ppm')
ax[1].set_xlim(100, 0)
ax[1].set_xlabel('Chemical Shift (ppm)')
ax[2].plot(output.spe['hz'], output.spe['real'])
ax[2].set_xlim(7500, 0)
ax[2].set_title('Spectrum in Hz')
ax[2].set_xlabel('Frequency (Hz)')
plt.show()
5. Advanced: Cross Polarization (CPMAS)¶
In a cross-polarization (CP) experiment, magnetization is typically transfered from an abundant nucleus (I) to a dilute via through-space (dipolar) coupling
Dipolar coupling and internuclear distance¶
The dipolar coupling constant \(b_{IS}\) between two spins \(I\) and \(S\) depends on their distance \(r\):
It falls off as \(1/r^3\), so nearby atoms have a much stronger coupling than distant ones.
In the spin system below, one 1H is coupled to two 13C nuclei with very different coupling strengths:
| Pair | \(b_{IS}\) (Hz) | Relative distance |
|---|---|---|
| H–C2 | −20 000 | closer |
| H–C3 | −5 000 | further |
The role of contact time¶
During the CP spin-lock magnetization flows from 1H to 13C at a rate governed by \(b_{IS}\), and thus a short contact time only allows magnetization to reach the strongly coupled (nearby) carbon, while a longer contact time gives enough time for the magnetization to travel further, also polarizing the weakly coupled (distant) carbon.
from simpyson.templates import CPMAS
# Spin system: one 1H coupled to two 13C at different distances.
# C2 has a strong dipolar coupling (-20000 Hz, closer to H).
# C3 has a weak dipolar coupling (-5000 Hz, further from H).
cp_spinsys = """
channels 1H 13C
nuclei 1H 13C 13C
shift 1 2.0p 0 0 0 0 0
shift 2 10p 0 0 0 0 0
shift 3 15p 0 0 0 0 0
dipole 1 2 -20000 0 0 0
dipole 1 3 -5000 0 0 0
"""
cp_calc = SimpCalc(
spinsys=cp_spinsys,
proton_frequency='400e6',
sw=20000,
spin_rate=10000,
start_operator="I1x",
detect_operator="I2p+I3p",
method="direct",
crystal_file="rep100",
gamma_angles=10,
verbose=0,
np=4096,
pulse_sequence='cp_mas',
lb=20,
)
# Short contact time: 500 μs — only the strongly coupled C2 (nearby) is polarized
cp_calc.pulse_sequence.update_parameters(pcp=500)
sim_short = cp_calc.run()
# Long contact time: 5000 μs — magnetization reaches the weakly coupled C3 (further away)
cp_calc.pulse_sequence.update_parameters(pcp=5000)
sim_long = cp_calc.run()
fig, ax = plt.subplots(figsize=(6, 4), constrained_layout=True)
ax.plot(sim_short.spe['ppm'], sim_short.spe['real'] / sim_short.spe['real'].max(), label='Short contact time (500 μs)')
ax.plot(sim_long.spe['ppm'], sim_long.spe['real'] / sim_long.spe['real'].max()+ 1, label='Long contact time (5000 μs)')
ax.set_xlabel('Chemical Shift (ppm)')
ax.set_ylabel('Normalised Intensity')
ax.set_xlim(20, 5)
ax.set_yticks([])
ax.legend()
ax.set_title('Effect of CP contact time on $^{13}C$ spectrum')
# Annotate peaks
ax.axvline(10, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
ax.axvline(15, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
ax.text(8.6, 0.15, 'C$_2$\n(10 ppm, close)', ha='center', fontsize=8, color='gray')
ax.text(14, 0.15, 'C$_3$\n(15 ppm, distant)', ha='center', fontsize=8, color='gray')
plt.show()
