Simulating slow-motion cw EPR spectra

This user guide explains how to simulate slow-motion cw EPR spectra of using EasySpin's function chili. It is assumed that you are familiar with the basics of MATLAB, esp. with structures.

This user guide contains the following topics: There are the following advanced topics:
Running the simulation

Slow-motion cw EPR spectra of S=1/2 systems are computed by the EasySpin function chili.

```chili(Sys,Exp)
```

It is called with two arguments. 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.

If no output argument is given, `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);
[Field,Spec] = chili(Sys,Exp);
```

The outputs `Field` and `Spec` are vectors containing the magnetic field axis 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 is

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

The first line defines the spin system, a nitroxide radical with anisotropic g and A tensors. The second line gives the rotational correlation time of the radical. The third line specifies experimental parameters, here only the microwave frequency (The magnetic field range is chosen automatically). The fourth line calls the simulation function, which plots the result. Copy and paste the code above to your MATLAB command window to see the graph.

Of course, the names of the input and output variables don't have to be `Sys`, `Exp`, `Field` 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 convenience, thoughout this tutorial, we will use short names like `Sys` and `Exp`.

Specifying the static parameters

The first input argument specifies the static parameters of the paramagnetic molecule. It is a MATLAB 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.

The field `Sys.Nucs` contains a string 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 multinuclear systems about details of the slow-motion simulation in that case.

`Sys.A` gives the hyperfine couplings in MHz (Megahertz), with three principal values on 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 13C 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 couplings are given in G (Gauss) or mT (Milltesla), so you have to convert these values. 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 mt2mhz:

```A_MHz = mt2mhz(A_mT);    % mT -> MHz conversion
A_MHz = mt2mhz(A_G/10);  %  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 3-element array (see the function erot).

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

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

Dynamic parameters

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

There are several alternative ways to specify the rate of isotropic rotational diffusion: either by specifying `tcorr`, the rotational correlation time in seconds

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

or by givings its base-10 logarithm

```Sys.logtcorr = -7;   % 10^-7 s = 100 ns
```

or by specifying the principal value of the rotational diffusion tensor (in s-1)

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

or by givings its base-10 logarithm

```Sys.logDiff = 9;   % 1e9 s^-1 = 1 ns^-1
```
`Diff` and `tcorr` are related by
```Diff = 1/6/tcorr;
```
The input field `Diff` can be used to specify an axial rotational diffusion tensor, by giving a 2-element vector with first the perpendicular and the the parallel principal value:
```Sys.Diff = [1 2]*1e8;  % in hertz
```

The `lw` field has the same meaning As for the other simulation functions garlic and pepper, the field `Sys.lw` can be used to specify a Gaussian and a Lorentzian broadening (FWHM, in mT)

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

For the reliability of the simulation algorithm we recommend to always use a small residual Lorentzian line width in `Sys.lw`

`chili` is also capable of simulating spectra including Heisenberg spin exchange. The effective exchange frequency (in MHz) is specified in `Sys.Exchange`, e.g.
```Sys.Exchange = 100;     % 100 MHz
```

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

```Nx = struct('g',[2.0088, 2.0061, 2.0027],'Nucs','14N','A',[16 16 86]);
Nx.tcorr = 1e-7; Nx.lw = [0 0.1]; Nx.Exchange = 100;
Exp = struct('mwFreq',9.5,'CenterSweep',[338, 10]);
chili(Nx,Exp);
```
The orienting potential

`chili` can also include an orienting potential in the simulation. It is specified in the field `lambda` in the spin system structure. Up to five coefficients can be given, λ2,0, λ2,2, λ4,0, λ4,2, and λ4,4, in that order. For example,

```Nx.lambda = [0.3, -0.2];
```

indicates λ2,0 = 0.3 and λ2,2 = -0.2.

Basic experimental settings

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
```

On many cw EPR spectrometers, the field range is specified using center field and sweep width, so `Exp.CenterSweep` is often the more natural choice.

`Exp.CenterSweep` and `Exp.Range` are only optional. If both are omitted, EasySpin tries to determine a field range large enough to accomodate the full spectrum. This automatic ranging works for most common systems, but fails in some 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 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, direct detection
Exp.Harmonic = 1; % first harmonic (default)
Exp.Harmonic = 2; % second harmonic
```

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.

More experimental settings

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. EasySpin can simulate both types of spectra:

```Exp.Mode = 'perpendicular'; % perpendicular mode (default)
Exp.Mode = 'parallel';      % parallel mode
```

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
```
Powder vs single orientation

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 ("director") is historically called the MOMD (microscopic order macroscopic disorder) model. This is equivalent to a powder spectrum. In `chili`, the powder spectrum is simulated whenever you specify an ordering potential. To get a single-crystal spectra, i.e. the slow-motion spectrum for the nitroxide attached to a rigid protein with a single orientation, give the crystal orientation in `Exp.CrystalOrientation`. `Exp.CrystalOrientaiton` contains the Euler tilt angle (in radians), betwen the lab frame (which is lab-fixed) and the frame of the orienting potential (which is fixed to the protein).

When `chili` performs a powder simulation, it takes the number of orientations to include from `Opt.nKnots`. Often, `Opt.nKnots` does not have to be changed from its default setting, but if the spectrum does not appear to be smooth, `Opt.nKnots` can be increased. Of course, this also increases the simulation time. For quick simulations, `Opt.nKnots` should be minimized.

Simulation parameters

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

The field `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 giving

```Opt.Verbosity = 1;     % log information
```

Another important option is `LLKM`. It specifies the number of orientational basis functions used in the simulation. For spectra in the fast and intermediate motion regime, the preset values don't have to be changed. However, close to the rigid limit, the default settings of `LLKM` might be too small. In that case, `LLKM` has to be increased, e.g.

```Opt.LLKM  = [24 20 10 10];
```

To see the values of `LLKM` that `chili` is using, set `Opt.Verbosity=1`, as described above.

Systems with more than one nucleus

`chili` is not capable of simulating a slow-motional cw EPR spectrum of systems with more than one nucleus by solving the Stochastic Liouville equation.

However, if the hyperfine coupling of one nucleus is significantly larger than those of the other nuclei, `chili` uses an approximate procedure: The slow-motional spectrum simulated using only the electron spin and the nucleus with the dominant hyperfine coupling is convoluted with the isotropic splitting pattern due to all other nuclei. This post-convolution technique gives resonable results.

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];
CuPc.logtcorr = -7.35;
Exp = struct('mwFreq',9.878,'CenterSweep',[330 120],'nPoints',5e3);
Opt.LLKM = [36 30 8 8];
chili(CuPc,Exp,Opt);
```
Frequency sweeps

`chili`, like the other cw EPR simulation functions `pepper` and `garlic`, does field sweeps 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 frequency range (in GHz) in `Exp.Range` or `Exp.CenterSweep`. You can also omit these, in which case `pepper` will determine an adequate range automatically.
• If you use `Sys.lw` or `Sys.lwpp`, make sure they are in MHz units. For a frequency sweep, these convolutional linewidth 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];
Nx.tcorr = 3e-9;
Exp.Field = 340;         % static field, in mT
Exp.Range = [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 ordering potential, you can run powder spectra, and you can adjust the basis size.