saffron(Sys,Exp)
saffron(Sys,Exp,Opt)
y = saffron(...)
[x,y] = saffron(...)
[x,y,info] = saffron(...)

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

Longitudinal relaxation time in microseconds. For simulations using fast the relaxation time has to be a scalar. Simulations that use the thyme method, it can be either a scalar or a matrix. If it is a scalar, the relaxation time is applied to all transitions. Transition specific relaxation times can be provided in the form of a matrix. See section Relaxation Times in the documentation on spidyan for more details. Relaxation times for transitions that are not defined (or 0) are automatically set to 10¹⁰ microseconds.

Transverse relaxation time in microseconds. For simulations using fast the relaxation time has to be a scalar. Simulations that use the thyme method, it can be either a scalar or a matrix. If it is a scalar, the relaxation time is applied to all transitions. Transition specific relaxation times can be provided in the form of a matrix. See section Relaxation Times in the documentation on spidyan for more details. Relaxation times for transitions that are not defined (or 0) are automatically set to 10¹⁰ microseconds.

Magnetic field (in mT) at which the experiment is performed.

Specifies a predefined pulse experiment or a user-defined pulse experiment.

Possible values for predefined experiments are '2pESEEM', '3pESEEM', '4pESEEM', 'HYSCORE', 'MimsENDOR'. See the page on predefined experiments for details.

User-defined sequences have to be given as cell array, with alternating entries for pulses (via a structure containing the pulse parameters) and numbers that give the inter-pulse delays in microseconds.

Pulse1.Type = 'rectangular';  % Pulse shape
Pulse1.tp = 0.032;            % Pulse length, microseconds
Pulse1.Flip = pi/2;           % Pulse flip angle, radians

Pulse2.Type = 'rectangular'; 
Pulse2.tp = 0.032;
Pulse2.Flip = pi;

Exp.Sequence = {Pulse1 0.2 Pulse2 0.1}; 

Ideal pulses (with zero length and therefore infinite excitation bandwidth) are defined by providing only a flip angle:

% Zero-length pulses
Pulse1.Flip = pi/2;
Pulse2.Flip = pi;

Exp.Sequence = {Pulse1 0.2 Pulse2 0.4};

The inter-pulse delays are defined as going from the end of one pulse to the beginning of the next.

EPR spectrometer frequency in GHz. For the the fast method,this needs only to be given if orientation selection is wanted (see Exp.ExciteWidth).

For the thyme method, this is required. All frequencies in the pulse definition (Pulse.Frequency) are defined relative Exp.mwFreq>:

Exp.mwFreq = 33.5;

Pulse.Frequency = [-250 250];  % corresponds to 33.25 - 33.75 GHz

A vector that contains the number of points in each dimension, e. g. [10 150] corresponds to 10 points in the first and 150 points in the second dimension.

Exp.nPoints = [10 150]; 	% 10 points along the first indirect dimension, 150 along the second

If omitted for predefined experiments, default values are assumed, see predefined experiments.

Dim1, Dim2,... provide the fields that are to be changed along the indirect dimensions (field nPoints). The first data point always uses the values defined initially in the experiment definition. All fields that appear in the pulse definition can be changed, e.g:

Exp.Dim1 = {'p2.Flip' pi/8};        % increments the flip angle of the second pulse by pi/8 each step

For free evolution events only the length can be changed:

Exp.Dim1 = {'d3' -0.1};             % decrements the length of the third delay by 100 ns each step

Several parameters can be changed in one dimension:

Exp.Dim1 = {'p2.Flip' pi/8; 'd3' -0.1}; 	% changes flip angle of pulse and duration of free evolution
Exp.Dim1 = {'p2.Flip,p3.Flip' pi/8};  		% flip angles of 2nd and 3rd pulse are simultaneously stepped  

For experiments that involve one or several moving pulses, the identifier Position can be used. This is only possible for pulses that are not the first or last event in the sequence. Pulses are allowed to cross, but must not overlap.

Exp.Dim2 = {'p2.Position' 0.1};   % moves the second pulse 100 ns back each step in the 2nd dimension

A list of increments can be used by providing vectors with the precomputed increments. All increments are always applied to the initial value of the field (as defined in the Exp) and are not related to each other. Hence, if the value in the experiment definition is desired as first data point of the indirect dimension, the fist element has to be zero:

Exp.nPoints = 4;

% this
Exp.Dim1 = {'d2' [0 0.1 0.2 0.3]};
% is equal to: 
Exp.Dim1 = {'d2' 0.1}; 

Complete freedom is given when it comes to providing the list of increments.

Exp.nPoints = 5;
Exp.Dim1 = {'d2' [0 0.1 0.3 0.65 -0.2]};   % increment of the second delay

However, the program checks that changing lengths of events do not lead to overlapping pulses at any acquisition point.

A special case is the field Par.Frequency in the pulse definition of a frequency-swept or of a monochromatic pulse that is defined by identical initial and final frequency. If one value is provided in Exp.DimX, this is going to be added to both elements in the Par.Frequency field - the pulse is moved in the frequency domain:

Exp.Dim1 = {'p1.Frequency' 5};   % changes both elements in Par.Frequency by 5 MHz per step

An equivalent input as the above is the following, where each the increment for each field in Par.Frequency is given:

Exp.Dim1 = {'p1.Frequency' [5 5]};   % changes both elements in Par.Frequency by 5 MHz per step

This syntax can also be used to change the sweep width:

Exp.Dim1 = {'p1.Frequency' [-5 5]};   			% increases sweep width by 10 MHz each step
Exp.Dim1 = {'p1.Frequency' [0 5]};    			% increment only the final frequency of the frequency sweep

A list of increments can be provided as well, by using ';' as separator:

Exp.Dim1 = {'p1.Frequency' [0 0; 10 10; 30 30; 80 80]};  % vector increment for shifting the frequency range of the pulse
Exp.Dim1 = {'p1.Frequency' [0 0; -10 0; 0 30; -80 80]};  % vector increment changing the excitation range of the pulse

For other pulse parameters that are defined by a vector (e.g., the order of an asymmetric hyperbolic secant (HS) pulse Pulse.n or the list of relatives amplitudes of a Gaussian cascade Pulse.A0), selected elements can be incremented by adding an index to the field name:

Exp.Dim1 = {'p1.A0(3)' 0.1; 'p1.A0(4)' -0.1};   % changes relative amplitudes of the third and fourth pulse in a Gaussian cascade
Exp.Dim1 = {'p1.n(2)' 2};    % increases order of the falling flank of a hyperbolic secant pulse by 2 each step

% Also possible for 'Frequency'
Exp.Dim1 = {'p1.Frequency(2)' 5};    % changes only the final frequency of the pulse
Exp.Dim1 = {'p1.Frequency' [0 5]};   % identical to the above and not a list of increments!

Adding an indirect dimension also allows for simultaneous phase cycling of two or more pulses (something that can not be achieved through the Exp.PhaseCycle structure):

Exp.nPoints = 4;
Exp.Dim1 = {'p2.Phase,p3.Phase' pi/4};   % changes the phase of the 2nd and 3rd pulse by pi/4 each step

In the above example, the output of spidyan will contain the individual transients from each phase cycling step and manual merging of the dimensions (with proper detection phase/sign) is required to obtain the phase-cycled signal.

This is an example for a two-dimensional experiment (e.g. for a two-pulse echo) to optimize two pulse parameters:

Exp.nPoints = [20 15];
Exp.Dim1 = {'p1.Frequency,p2.Frequency' [-2.5 2.5]}; % increase excitation band by 5 MHz each step of both pulses
Exp.Dim1 = {'p1.tp' 2; 'p2.tp' 1}; % step pulse lengths, in µs

An Nx3 array that specifies the sample orientations for which the EPR spectrum should be computed. Each row of SampleFrame contains the three Euler rotation angles that transform the lab frame to the sample/crystal frame.

Exp.SampleFrame = [0 0 0];                   % sample/crystal frame aligned with lab frame
Exp.SampleFrame = [0 pi/2 0];                % sample/crystal frame tilted relative to lab frame
Exp.SampleFrame = [0 pi/2 pi/4];             % sample/crystal frame tilted relative to lab frame
Exp.SampleFrame = [0 0 0; 0 pi/2 pi/4];      % two samples/crystals

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

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.

Cell array that specifies the phase cycle, with one entry for each pulse. To cycle the phase of a pulse, provide an array with each phase step in a separate row. In each row the first element is the pulse phase angle (in radians), and the second element the detection phase (+1 or -1 or 1i or -1i).

PC2 = [0, +1; pi, -1];   % step 1: pulse phase angle = 0, detection phase = +1
                         % step 2: pulse phase angle = pi, detection phase = -1

Exp.PhaseCycle = {[0 +1] PC2};  % phase cycle only second pulse; use [0 1] for non-cycled pulses
Exp.PhaseCycle = {[] PC2};      % alternative: [] for non-cycled pulses

The phase cycle is built from the individual phase steps for each pulse in a nested fashion, where each pulse is cycled separately. To step more than one pulse simultaneously, use an additional indirect phase dimension to the experiment (see Dim).

Optional, can be 'linear' or 'circular', default is 'linear'. If mwPolarization is 'circular', the excitation operator also uses the imaginary part of the waveform. The default excitation then becomes real(IQ) Ŝ_x + imag(IQ) Ŝ_y.

Exp.mwPolarization = 'circular';     % mwPolarziation set to circular for all pulses

Keep in mind, that using 'circular' does not allow saffron to use some of its speed up tricks and can lead to significantly longer simulation times.

The field DetWindow controls detection and is given in microseconds: The zero-time of the detection is taken as the end of the last element in Exp.Sequence. It only works with the thyme method. The detected signal corresponds to ⟨Ŝ_-⟩, which means the real part corresponds to ⟨Ŝ_x⟩ and the imaginary part to ⟨-Ŝ_y⟩.

To do single point detection, it is sufficient to tell DetWindow the time of the detection point. E.g. to detect at the center of the echo after a two pulse experiment:

tau = 0.5
Exp.Sequence = {pulse90 tau pulse180 tau};   	% a two-pulse echo

Exp.DetWindow = 0;	% detect a single point at the end of the second delay (center of the echo)

The detection position can be set to anywhere at the end of the sequence, as long as it does not overlap with a pulse:

Exp.DetWindow =  0.2;	% detect a single point 200 ns after the second delay
Exp.DetWindow = -0.2;	% detect a single point 200 ns before the end of the second delay

In order to detect a transient, DetWindow has to be provided with a start and end value:

tau = 0.5
Exp.Sequence = {pulse90 tau pulse180 tau};   	% a two-pulse echo

Exp.DetWindow = [-0.1 0.1];	% detect a from -100 ns to 100 ns centered around of the end of tau

Of course, it is also possible to detect a transient anywhere after the last pulse, as long is it does not overlap with it:

Exp.DetWindow = [0 0.3];	% start detection at the end of the last event and detect for 300 ns
Exp.DetWindow = [-0.3 -0.1];	% start detection 300 ns before the end of the last event and detect for 200 ns

If you want to detect an FID or have a fixed experiment length (e.g. 3p ESEEM) you can also write:

Exp.Sequence = {pulse90 0.5 pulse90 0.8 pulse90};   	% a three-pulse ESEEM sequence (the last element is a pulse)

Exp.DetWindow = 0.5;	% detect a single point after 500ns after the last pi/2 pulse (center of echo)

Exp.DetWindow = [0.250 0.750] % detect a transient, centered around the echo

Exp.DetWindow = [0 0.100]    % detect an FID after the last pulse

The signal is detected in the simulation frame and hence contains an oscillating component. By translating/shifting the frequency it is possible to center the signal around 0 (baseband) or any other desired frequency. The frequency shift has to given as absolute frequency :

Exp.DetFreq = 33.5; % shifts detected signal by -33.5 GHz

If no DetFreq is provided, Exp.mwFreq is taken as frequency for down-conversion.

Phase of detection in radians, optional. The default value is 0.

Exp.ResonatorFrequency = 33.5; 	% resonator center frequency in GHz
Exp.ResonatorQL = 300;	  	% loaded Q-value

Exp.FrequencyResponse = [Frequency; TransferFunction] 	% transfer function and its correspond frequency axis

Center frequency of resonator in GHz.

The loaded Q-value.

This gives the frequency response of the resonator, in the form Par.FrequencyResponse = [Frequency; TransferFunction] with the (possibly experimental) resonator transfer function in TransferFunction and the corresponding frequency axis in Frequency (in GHz). A complex transfer function input in Par.FrequencyResponse is used directly in the bandwidth compensation. A real transfer function input is assumed to correspond to the magnitude response, and the associated phase response is estimated.

If FrequencyResponse is given, the alternative input fields ResonatorFrequency and ResonatorQL are ignored.

Optional, can be 'simulate' or 'compensate'. By default the effect of the resonator on the signal is simulated. If set to 'compensate' the pulse shape is adapted such that it compensates for the resonator profile by using a uniform critical adiabaticity criterion.

A number giving the number of orientations (knots) for which spectra are explicitly calculated. GridSize 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.

Optional, allows to switch simulation mode. Can be 'fast' for the fast algorithm or 'thyme' the general approach. By default, saffron choses 'fast' if possible.

Determines whether to return the total spectrum or a list of subspectra. If set to '', the total spectrum is returned. If set to 'components', spec is a matrix with the subspectra of all components (including isotopologues). Each row in spec is one subspectrum.

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:

• EndorMethod = 0: This is a sum-over-transitions method that emulates rf pulses by swapping populations of adjacent nuclear spin levels. For one nucleus, it is roughly equivalent to EndorMethod=1, but faster. For multiple nuclei, it gives wrong results: With ProductRule=1, cross suppression effects are not modeled, and with ProductRule=0, peak positions are wrong.
• EndorMethod = 1: This is a sum-over-transitions method that applies bandwidth-limited rf pi pulses (using single-transition operators in the eigenbasis of the nuclear sub-Hamiltonians) on each nuclear transition in turn. It is able to reproduce inter-nuclear cross suppression effects (implicit triple). All rf pulses are modeled as 180 degree pulses on all allowed nuclear transitions, independent of the nature of the nucleus.
• EndorMethod = 2: This is an alternative method. It uses a very slow brute-force rf sweep approach: It loops over every frequency on the rf axis and calculates the echo amplitude using the same rf pulse operators as EndorMethod=1.

The default method is 0 for a single nucleus or multiple nuclei with product rule, otherwise it is 1.

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.

Optional, frequency of the (rotating) simulation frame in GHz. If SimFreq is not provided, a suitable simulation frame is automatically selected. No matter the choice of the simulation frame, all other frequencies must still be given in the lab frame - the program handles all required frequency shifts. For example, for running a Q-band simulation without using the default simulation frame, but one that rotates at 30 GHz:

Exp.mwFreq = 34.4;	        % Set carrier frequency to 34.4 GHz
Pulse.Frequency = [-50 50]; % Pulse with bandwidth of 100 MHz, sweeps from 34.35 to 34.45 GHz 

Opt.SimFreq = 30;  	      % Changes into a rotating frame at 30 GHz

To force a simulation in the labframe SimFreq has to be set to zero:

Opt.SimFreq = 0;  	      % Propagation in the labframe

Changing the simulation frame allows for using a larger time step, which in turn can strongly reduce computation time. Keep in mind, that if a time step is provided with Opt.IntTimeStep it is not adapted automatically. To have the program adapt the time step to the simulation frequency, remove the field Opt.IntTimeStep .

Integration time step in microseconds, optional. If no time step is provided, the program computes a suitable time step (taking into account the highest frequency that is being used for the pulse definition Exp.Frequency). Accurate results are obtained for an integration time step that corresponds to about 50 data points per oscillation or 1/50th of the largest time step size to fulfil the Nyquist criterion.

User provided integration time steps must at least fulfil the Nyquist criterion, or an error is thrown. For integration time steps larger than the recommended ones, a warning is printed to the command line.