This chapter explains how to simulate slow-motion cw EPR spectra. "Slow-motion" indicates that each spin center (radical or transition metal complex) in the sample does not have a fixed orientation, but reorients/tumbles randomly. This motion is called rotational diffusion. Depending on the rate of rotational diffusion, different simulation functions need to be used: If the rotational motion is frozen, the sample is in the rigid limit, and pepper is the simulation function to use. If the rate of rotational diffusion is very fast, then garlic should be used. For intermediate tumbling rates, the regime is called the "slow-motion regime". For this regime, EasySpin's function chili should be used. Here, we show how to use chili.

• Running the simulation
• Static parameters
• Dynamic parameters
• The orienting potential
• Basic experimental settings
• More experimental settings
• Orientational averaging
• Simulation parameters
• Spin centers with many nuclei
• Frequency-swept spectra

Slow-motion cw EPR spectra of S=1/2 systems are computed by the EasySpin function chili. Just as the other simulation functions (pepper, garlic, etc.), it is called with two or three arguments:

chili(Sys,Exp)
chili(Sys,Exp,Opt)

The first argument Sys tells chili all about the static and dynamic parameters of the spin system. The second argument Exp gives the experimental parameters. The third argument Opt gives calculation options.

Without output arguments, chili plots the computed spectrum. But it can also return one or two outputs. (Don't forget the semicolon at the end of the line to suppress output to the command window.)

Spec = chili(Sys,Exp);
[x,Spec] = chili(Sys,Exp);

The outputs x and Spec are vectors containing the magnetic field axis (in millitesla) and the spectrum, respectively. If these are requested, chili does not plot the spectrum.

Doing a simulation only requires a few lines of code. A simple example of a tumbling nitroxide radical is

clear
Sys.g = [2.008,2.006,2.003];
Sys.Nucs = '14N';
Sys.A = [20,20,85];    % MHz
Sys.tcorr = 4e-8;      % seconds
Exp.mwFreq = 9.5;      % GHz
chili(Sys,Exp);

The first few lines define the spin system, a nitroxide radical with anisotropic g and A tensors. Then the rotational correlation time of the radical is set. Finally, experimental parameters are defined, here just the microwave frequency (the magnetic field range is chosen automatically). The last line calls the simulation function, which plots the result. Copy and paste the code above to your MATLAB command window to see the result.

Of course, the names of the input and output variables don't have to be Sys, Exp, x, and Spec. You can give them any name you like, as long as it is a valid MATLAB variable name, e.g., FremySaltSolution or QbandExperiment. For economy, throughout this tutorial, we will use short names like Sys and Exp.

The first input argument to chili specifies the static parameters of the paramagnetic molecule. It is a structure with various fields giving values for the spin system parameters.

Sys.g = [2.008,2.006,2.003];
Sys.Nucs = '14N';
Sys.A = [20,20,80];  % MHz

The first line defines the g tensor of the spin system via its three principal values. chili always assumes a single unpaired electron spin S=1/2 by default if nothing is given about the spin S. For other values of S, Sys.S must be given. For example, if you want to simulate a triplet, use Sys.S=1.

The field Sys.Nucs contains a character array giving all the magnetic nuclei in the spin system, a nitrogen-14 in the above example. Use a comma-separated list of isotope labels to give more than one nucleus. E.g., Sys.Nucs = '14N,1H,1H' specifies one nitrogen and two different protons. See the section on multi-nuclear systems about details of the slow-motion simulation in that case.

Sys.A gives the hyperfine coupling tensors in MHz (megahertz), with up to three principal values in a row for each of the nuclei listed in Sys.Nucs. The following defines a hydrogen atom with a 10 MHz coupling to the unpaired electron and a ¹³C atom with a 12 MHz coupling.

Sys.Nucs = '1H,13C';
Sys.A = [10; 12]; % MHz

Remember that chili (and other EasySpin functions, too), take the hyperfine coupling values to be in MHz. Often, values for hyperfine coupling constants are given in G (gauss) or mT (millitesla), so you have to convert these values to MHz. For g = 2.00232, 1 G corresponds to 2.8025 MHz, and 1 mT corresponds to 28.025 MHz. The simplest way to convert coupling constants from magnetic-field units to MHz is to use the EasySpin function unitconvert:

A_MHz = unitconvert(A_mT,'mT->MHz');    % mT -> MHz conversion
A_MHz = unitconvert(A_G/10,'mT->MHz');  %  G -> MHz conversion (1 G = 0.1 mT)

The orientations of the tensors relative to the molecular frame are defined in terms of Euler angles in a 3-element array (see the function erot and the section on frames).

Sys.gFrame = [0 0 0];    % Euler angles for g tensor frame
Sys.AFrame = [0 pi/4 0]; % Euler angles for A tensor frame

For more details on these static parameters of the spin system, see the documentation on spin systems.

The spin system structure also collects parameters relating to the dynamics of the paramagnetic molecules. chili of Brownian diffusion

The most important parameter is the rate of rotational diffusion/tumbling. It can be isotropic or anisotropic (i.e. diffusion can be faster around some axes than around others). There are several alternative ways to specify the rate of isotropic rotational diffusion: either by specifying Sys.tcorr, the rotational correlation time in seconds

Sys.tcorr = 1e-7;   % 10^-7 s = 100 ns

by giving its base-10 logarithm

Sys.logtcorr = -7;   %  = log10(Sys.tcorr)   10^-7 s = 100 ns

or by specifying the principal value of the rotational diffusion tensor (in rad² s^-1, or equivalently s^-1)

Sys.Diff = 1e9;  % 1e9 rad^2 s^-1 = 1 rad^2 ns^-1

or by giving its base-10 logarithm

Sys.logDiff = 9;   % = log10(Sys.Diff)   1e9 rad^2 s^-1 = 1 rad^2 ns^-1

Diff and tcorr are related by

Diff = 1/6./tcorr;

Sys.Diff = [1 2]*1e8;  % in rad^2 s^-1

The Sys.lw field has the same meaning as for the other simulation functions garlic and pepper. It can be used to specify an additional Gaussian and/or Lorentzian broadening (FWHM, in mT):

Sys.lw = [0.05 0.01];   % [GaussianFWHM, LorentzianFWHM] in mT

These line broadenings are convolutional, i.e. they broaden all regions of the spectrum equally. For the reliability of the slow-motion simulation algorithm in chili, it is recommended to always use a small residual Lorentzian line width in Sys.lw.

Sys.Exchange = 100;     % microseconds^-1

A short example of a nitroxide radical EPR spectrum with exchange narrowing is

clear
Nx.g = [2.0088, 2.0061, 2.0027];
Nx.Nucs = '14N';
Nx.A = [16 16 86];  % MHz
Nx.tcorr = 1e-7;    % seconds
Nx.lw = [0 0.1];    % mT
Nx.Exchange = 30;   % microseconds^-1
Exp.mwFreq = 9.5;
Exp.CenterSweep = [338 10]; % mT
chili(Nx,Exp);

chili can also include an orientational potential-energy function in the simulation. This is used to model hindered rotational diffusion, i.e. when not all orientations of the spin center relative to its immediate environment (protein, membrane, etc) are equally accessible energetically.

The orientational potential is specified in the field Sys.Potential in the spin system structure. The potential is expressed as a linear combination of Wigner function, with the expansion coefficients λ_LMK and the function quantum numbers L, M, and K. For each expansion term, specify L, M, K, and λ_LMK in a row. Here are some examples:

Sys.Potential = [2 0 0 0.3];              % one term
Sys.Potential = [2 0 0 0.3; 2 0 2 -0.2];  % two terms

The orientational potential gives rise to an non-uniform orientational distribution of the spins. The simplest ordering term is λ_2,0,0. Its symmetry corresponds to a d_z2 orbital. If it is positive, it indicates preferential orientation of the molecular z axis along the external magnetic field.

When an ordering potential for the spin center is present, the simulation has to take into account whether the immediate environment responsible for the ordering (protein, membrane, etc) is ordered or disordered. In the former case, a single-orientation simulation is sufficient, whereas in the latter case, an orientational average over all possible orientations of the protein or membrane is needed - see section on orientational averaging.

The second input argument, Exp, collects all experimental settings. Just as the spin system, Exp is a structure containing several fields.

Microwave frequency. To simulate an EPR spectrum, EasySpin needs at a minimum the spectrometer frequency. Put it into Exp.mwFreq, in units of GHz.

Exp.mwFreq = 9.385;  % X-band
Exp.mwFreq = 34.9;   % Q-band

Field range. There are two ways to enter the magnetic-field sweep range. Either give the center field and the sweep width (in mT) in Exp.CenterSweep, or specify the lower and upper limit of the sweep range (again in mT) in Exp.Range.

Exp.CenterSweep = [340 80]; % in mT
Exp.Range = [300 380];      % in mT

Exp.CenterSweep and Exp.Range are only optional. If both are omitted, EasySpin tries to determine a field range large enough to accommodate the full spectrum. This automatic ranging works for many common systems, but fails in more complicated situations. EasySpin will issue an error when it fails.

Points. By default, pepper computes a 1024-point spectrum. However, you can change the number of points to a different value using

Exp.nPoints = 5001;

You can set any value, unlike some older EPR spectrometers, where often only powers of 2 are available (1024, 2048, 4096, 8192).

Harmonic. By default, EasySpin computes the first-harmonic absorption spectrum, i.e. the first derivative of the absorption spectrum. By changing Exp.Harmonic, you can request the absorption spectrum directly or the second-harmonic (second derivative) of it.

Exp.Harmonic = 0; % absorption spectrum
Exp.Harmonic = 1; % first harmonic (default)
Exp.Harmonic = 2; % second harmonic

Often, it is useful to look at the simulated absorption spectrum (Exp.Harmonic=0) first to develop an understanding of spectral features.

Modulation amplitude. If you want to include effects of field modulation like overmodulation, use Exp.ModAmp

Exp.ModAmp = 0.2; % 0.2 mT (2 G) modulation amplitude, peak-to-peak

Time constant. To include the effect of the time constant, apply the function rcfilt to the simulated spectrum. Although many existing spectrometers have this electronic filtering, it is considered outdated and should be avoided in experiments.

For more advanced spectral simulations, EasySpin offers more possibilities in the experimental parameter structure Exp.

Mode. Most cw EPR resonators operate in perpendicular mode, i.e., the oscillating magnetic field component of the microwave in the resonator is perpendicular to the static field. Some resonators can operate in parallel mode, where the microwave field is parallel to the static one. Although EasySpin can simulate both types of spectra in the rigid limit, it does not support parallel mode for slow-motional spectra.

Temperature. The polarizing effects of low sample temperatures can also be included in the simulation by specifying the temperature:

Exp.Temperature = 4.2; % temperature in kelvin

With this setting, EasySpin will include the relevant polarization factors resulting from a thermal equilibrium population of the energy levels. For S=1/2 systems, it is not necessary to include the temperature. However, it is important in high-spin systems with large zero-field splittings, and in coupled spin systems with exchange couplings.

Microwave phase. Occasionally, the EPR absorption signal has a small admixture of the dispersion signal. This happens for example when the microwave phase in the reference arm is not absolutely correctly adjusted. EasySpin can mix dispersion with absorption if a Lorentzian broadening is given:

Sys.lwpp = [0.2 0.01];           % Lorentzian broadening (2nd number) required

Exp.mwPhase = 0;                 % pure absorption
Exp.mwPhase = pi/2;              % pure dispersion
Exp.mwPhase = 3*pi/180;          % 3 degrees dispersion admixed to absorption

In a frozen solution of spin-labelled protein, the protein molecules assume all possible orientations. For slow-motion spectra, this orientational distribution has to be taken into account if a orienting potential is present. If not, it is sufficient to compute only one orientation, as the spectra from all orientations are identical.

The summation of slow-motion spectra over all possible orientations of an immobile protein or membrane ("director") is historically called the MOMD (microscopic order macroscopic disorder) model. This is equivalent to a powder average. In chili, the orientational average is calculated whenever you specify an orientational potential. To get a single-orientation spectra, i.e. the slow-motion spectrum for the nitroxide attached to a rigid protein with a single protein orientation, use the field Exp.SampleFrame to specify the sample orientation via the Euler tilt angle (in radians) between the lab frame (which is lab-fixed) and the frame of the orienting potential (which is fixed to the protein).

When chili performs an orientational average, it takes the number of orientations to include from Opt.GridSize. Often, Opt.GridSize does not have to be changed from its default setting, but if the spectrum does not appear to be smooth, Opt.GridSize can be increased. Larger values of Opt.GridSize increases the simulation time. It might also be necessary to change the grid symmetry in Opt.GridSymmetry (to e.g. 'D2h' or 'Ci').

chili is not able to automatically determine values of Opt.GridSize and Opt.GridSymmetry that yield a converged spectrum. It is therefore important to increase Opt.GridSize until the spectrum doesn't change anymore, and to initially use the most general symmetry setting Opt.GridSymmetry = 'Ci' to assure a correct and converged spectrum.

The third input structure, Opt, collects parameters related to the simulation algorithm.

The most important field in this structure is Opt.LLMK. It contains 4 numbers that determine the number N of orientational basis functions used in the simulation - the larger the numbers in LLMK the larger is N. N determines how accurately the rotational motion is modelled: The larger N, the more accurate the simulation. Typically, the spectrum converges as a function of N. chili uses preset values in Opt.LLMK which typically work well for X-band nitroxide spectra in the fast and intermediate motion regime without ordering. However, in the presence of ordering, for slow motions, for other spin systems, and for other spectrometer frequencies, the default settings of Opt.LLMK are too small. In that case, the four values in Opt.LLMK have to be increased, e.g.

Opt.LLMK  = [24 20 10 10];  % moderately high values

All numbers must be positive or zero. The first number is the maximum even L, the second is the maximum odd L, and the third and fourth number are the maximum M and K, respectively. The latter two numbers cannot be larger than the maximum L.

One approach to optimize the basis size is to start from [2 0 0 0] and increase each of the four numbers in turn, run the simulation, and see if the spectrum changes. Once the spectrum doesn't change anymore, the basis size large enough. If your spin system or rotational parameters change, you need to test the basis size again for convergence.

To see the values of Opt.LLMK that chili is using, set Opt.Verbosity=1, as described next.

The field Opt.Verbosity specifies whether chili should print information about its computation into the MATLAB command window. By default, its value is set to 0, so that chili is silent. It can be switched on by setting it to a non-zero value.

Opt.Verbosity = 0;     % don't print any information
Opt.Verbosity = 1;     % print information
Opt.Verbosity = 2;     % print more information

The Stochastic Liouville equation solver in chili is capable of simulating slow-motional cw EPR spectra of systems with multiple nuclei. However, the computational cost increased rapidly with the number of spins, to a degree that simulations with more than three or four nuclei might be unbearably slow or may run out of memory.

There is a workaround for situations where the hyperfine coupling of one or two nuclei is significantly larger than those of all other nuclei. In that case, chili can be instructed via Opt.PostConvNuclei to use an approximate procedure: The slow-motional spectrum is simulated using only the electron spin and the nuclei with the dominant hyperfine couplings, and the resulting spectrum is convoluted with the isotropic splitting pattern due to all other nuclei. This post-convolution technique often gives reasonable results very close to the exact solution.

A simple example is

CuPc = struct('g',[2.0525 2.0525 2.1994],'Nucs','63Cu,14N','n',[1 4]);
CuPc.A = [-54 -54 -608; 52.4 41.2 41.8]; % MHz
CuPc.logtcorr = -8.3;
Exp = struct('mwFreq',9.878,'CenterSweep',[330 120],'nPoints',5e3);
Opt.LLMK = [36 30 8 8]; % fairly large basis
Opt.PostConvNucs = 2;  % indicate to use post-convolution for 14N
chili(CuPc,Exp,Opt);

chili, like the other cw EPR simulation functions pepper and garlic, simulates field-swept spectra by default. However, you can use it to simulate frequency-swept spectra as well.

For this, all you need to do is the following:

• Give a static magnetic field (in mT) in Exp.Field. Make sure you do not set Exp.mwFreq, otherwise EasySpin does not know what to do.
• Give a microwave frequency range (in GHz) in Exp.mwRange or Exp.mwCenterSweep. You can also omit these, in which case pepper will try to determine an adequate range automatically.
• If you use Sys.lw or Sys.lwpp, make sure they are in units of MHz. For a frequency sweep, these convolutional line width parameters are understood to be in MHz (and not in mT, as they are for field sweeps).

Here is an example of a frequency-swept slow-motion spectrum of a nitroxide radical:

clear
Nx.g = [2.008 2.006 2.002];
Nx.Nucs = '14N';
Nx.A = [20 20 100];        % MHz
Nx.tcorr = 1e-9;           % seconds
Nx.lw = [0 5];             % MHz!
Exp.Field = 340;           % static field, in mT
Exp.mwRange = [9.3 9.7];   % frequency range, in GHz
chili(Nx,Exp);

By default, chili returns the absorption spectrum (Exp.Harmonic=0) when you simulate a frequency-swept spectrum. To get the first or second derivative, change Exp.Harmonic to 1 or 2. Note however that Exp.ModAmp is not supported for frequency sweeps.

All other capabilities of chili apply equally to frequency sweep and to field sweeps. For example, you can simulate multi-component spectra, you can use an orientational potential, you can perform orientational averaging, and you can adjust the basis size.