Calculation of pulse EPR spectra.

y = saffron(...)
[x,y] = saffron(...)
[x,y,p] = saffron(...)
[x1,x2,y,p] = saffron(...)

See also the examples on how to use saffron.


This function calculates pulse EPR (ESEEM and ENDOR) spectra of powder and single crystals.

The output contain the abscissa data in x (time in microseconds or frequency in MHz) and the simulated data in y (time domain trace or ENDOR spectrum). For 2D experiments such as HYSCORE, the two axes are returned in x1 and x2. For ESEEM simulations, p contains the frequency abscissa (in MHz) in p.f and the spectrum obtained by Fourier transform from the time domain data in p.fd.

If you don't request any output, saffron plots the simulated data.

The three input arguments to the function are

Sys is a spin system structure. Fields available in Sys include all needed for the construction of the spin Hamiltonian. Line broadening parameters used by other simulation functions (lw, lwpp, gStrain, etc.) are not recognized, except HStrain. HStrain is used in excitation window computations (see Exp.ExciteWidth) when orientation selection is wanted.

saffron supports any spin system with one electron spin (arbitrary S) and any number of nuclei.

If no orientation selection is required, then even the g tensor (and the microwave frequency) can be omitted. Only the nuclear parameters (and the field) need to be given:

Sys.Nucs = '14N';
Sys.A_ = [0.2 0.3];
Sys.Q = [-1 -1 2]*0.1;

You can provide the transverse and longitudinal relaxation times in the spin system structure:

An array [T1 T2] with two numbers, the longitudinal relaxation time constant T1 and the transverse relaxation time constant T2.

Exp contains the experimental parameters, most importantly the magnetic field and the pulse sequence.

Magnetic field (in mT) at which the experiment is performed.
Specifies the pulse sequence, '2pESEEM', '3pESEEM', '4pESEEM', 'HYSCORE', or 'MimsENDOR'.
Number of points in the ESEEM time trace or in the ENDOR spectrum. If only one number is given and a two-dimensional experiment is simulated, the number is used in both dimensions. If you give two numbers, each refers to one dimension. E.g. [31 256] indicates 31 points in the first dimension and 256 points in the second dimension. If not given, default values are used.
For ESEEM, time increment (dwell time), in microseconds. For 1D experiments, this should be a single number. For 2D sequences, dt applies to both dimensions if it's only one number, alternatively one number for each dimension can be given, e.g. dt = [0.024 0.008].
For three-pulse ESEEM, four-pulse ESEEM, HYSCORE and Mims ENDOR, this gives the value of τ, in microseconds. For two-pulse ESEEM, it is the starting τ value.
For three-pulse ESEEM and four-pulse ESEEM, this gives the starting T value, in microseconds. For the other sequences, this value has no effect.
For HYSCORE, these give the starting values of t1 and t2, in microseconds.
Contains lower and upper limit of the frequency range for ENDOR. Values should be in MHz, e.g. Exp.Range=[0 30].

For orientation selection, the following additional parameters are needed.

EPR spectrometer frequency in GHz. Needs only to be given if orientation selection is wanted.
The microwave excitation bandwidth in MHz (responsible for orientation selection). The excitation profile is assumed to be Gaussian, and ExciteWidth is its FWHM. The default is infinity. To obtain the full excitation with for a given orientation, ExciteWidth is combined with HStrain from the spin system structure.

For user-defined pulse experiments, the following fields are

List of pulse flip angles in multiples of 90°, e.g. [1 2] for two-pulse ESEEM and [1 1 2 1] for HYSCORE.
List of pulse lengths, in microseconds. If not given, all pulses are assumed to be infinitely short (ideal). If given, saffron integrates the signal over a small offset distribution. See Opt.nOffsets and Opt.lwOffset. Some of the values in tp can be zero, in which case ideal pulses are used. E.g. [0.200 0 0] is a three-pulse sequence with one selective finite-length pulse followed by two ideal pulses.
List of pulse phases, in multiples of 90°. Specifically, 0 is x, 1 is y, 2 is -x, 3 is -y.
List of initial delays, in microseconds. E.g. [0.2 0.2] for two-pulse ESEEM with starting τ of 200 ns, or similarly [0.2 0.1 0.1 0.2] for HYSCORE. The delays are defined to go from the end of one pulse to the beginning of the next (unlike in Bruker spectrometers).
Incrementation information for every inter-pulse delay. 0 indicates no incrementation, 1 increment along first dimension, 2 increment along second dimensions. Example: [0 1 2 0] for HYSCORE and [1 1] for two-pulse ESEEM.
Coherence filter, one character for each inter-pulse delay. '+' stands for electron coherence order +1, '-' for -1, 'a' for 0 (alpha), 'b' for 0 (beta), '0' for 0 (alpha or beta), 1 for +1 or -1, and '.' for anything. Examples: '.ab.' selects the coherence transfer pathways in HYSCORE that leads to alpha/beta cross peaks.

To simulate single crystals, use

An Nx3 array that specifies single-crystal orientations for which the EPR spectrum should be computed. Each row of CrystalOrientation contains the three Euler rotation angles that transform the crystal frame (C) to the lab frame (L). If CrystalOrientation is empty or not specified, the full powder spectrum is computed.
Exp.CrystalOrientation = [0 0 0];              % single crystal, crystal z axis aligned with B0
Exp.CrystalOrientation = [0 pi/2 0];           % single crystal, crystal z axis perpendicular to B0
Exp.CrystalOrientation = [0 0 0; 0 pi/2 0];    % two crystals
Exp.CrystalOrientation = [];                   % powder
Specifies the symmetry of the crystal. You can give either the number of the space group (between 1 and 230), the symbol of the space group, or the symbol for the point group of the space group (in either Schönflies or Hermann-Mauguin notation).

Exp.CrystalSymmetry = 'P21/c'; % space group symbol
Exp.CrystalSymmetry = 11;      % space group number (between 1 and 230)
Exp.CrystalSymmetry = 'C2h';   % point group, Schönflies notation
Exp.CrystalSymmetry = '2/m';   % point group, Hermann-Mauguin notation

When CrystalSymmetry is given, all symmetry-related sites in the crystal are included in the calculation. If CrystalSymmetry is not given, space group 1 (P1, point group C1, one site per unit cell) is assumed.

The fields in the structure Opt specify parameters for the simulation algorithm.

A number giving the number of orientations (knots) for which spectra are explicitly calculated. nKnots gives the number of orientations on quarter of a meridian, i.e. between θ = 0 and θ = 90°. The default value is 31, corresponding to a 3° spacing between orientations. For highly anisotropic spectra, esp. for HYSCORE, the value often has to be increased to 181 (0.5° spacing) or beyond.
0 (default) or 1. Determines whether saffron generates the spectrum by binning all the peaks in the frequency domain or by evolution of all the complex exponentials in the time domain. The frequency-domain binning method is very fast and is therefore used as the default. However, it involves small rounding of peak positions, which can in some cases lead to imperfect phase interferences and small artifacts in the spectrum. The time-domain method is significantly slower, but accurate.
Expansion factor used in the simulation, should be between 0, 1, 2, 3 or 4. The higher, the more accurate is a simulation, but the slower it becomes, especially for 2D simulations. Default values are 4 for 1D and 2 for 2D.
Determines whether product rule is used or not (1 or 0). By default, it is not used, but simulations with spin systems with more than two nuclei it might run faster with the product rule. The spectral result is independent of this choice.
Determines which method to use to simulate ENDOR spectra. There are three methods:
Number of points for the frequency offset integration in the case of finite-length pulses. Typical values are between 10 and 100, but should be determined for each case individually.
FWHM (in MHz) of the Gaussian distribution of offset frequencies use in the offset integration in the case of finite-length pulses. Typically, this is about the inverse of the length of the first pulse in a pulse sequence, e.g. 100 MHz for a 10ns pulse.
A 1 indicates that the HYSCORE spectrum should be plotted with a logarithmic intensity scale. If 0 (the default), a linear scale is used.
Apodization window used before FFT. See apowin for details.
The factor by which the time domain signal array should be padded with zeros before FFT. E.g. with ZeroFillFactor=4, a 256-point array is padded to 1024 points.

For both ESEEM and ENDOR, saffron uses matrix-based methods similar to those employed by Mims (1972) to compute frequencies and amplitudes of all peaks. With these peaks, a spectrum histogram is constructed, from which the time-domain signal is obtained by inverse Fourier transform.

For the pre-defined sequences, saffron assumes ideal pulses with standard flip angles (two-pulse ESEEM: 90°-180°, three-pulse ESEEM 90°-90°-90°, HYSCORE 90°-90°-180°-90°).

For systems with several nuclei, saffron by default simulates without using product rules, but can employ them if wanted (see Options).

For high-electron spin systems, all terms in the zero-field splitting, even the nonsecular ones, are taken into account.

To generate the spectrum from the time-domain signal, saffron performs (1) baseline correction, (2) apodization with a Hamming window, (3) zero-filling, and (4) FFT.

All the theory is described in

See also

nucfrq2d, salt, pepper