pepper(Sys,Exp);
pepper(Sys,Exp,Opt);
spec = pepper(...);
[B,spec] = pepper(...);
[B,spec,info] = pepper(...);
[nu,spec] = pepper(...);
[nu,spec,info] = pepper(...);

Specifies the microwave excitation mode. Possible settings are

Exp.mwMode = 'perpendicular';  % default
Exp.mwMode = 'parallel';
Exp.mwMode = {k, pol};

Resonator experiments:
For conventional EPR experiments with linearly polarized microwave in a resonator, use 'perpendicular' (default) or 'parallel'. In the perpendicular mode, the microwave magnetic field B₁ is oscillating along the laboratory x axis (xL), perpendicular to the external static magnetic field B₀. In the parallel mode, it is oscillating along the laboratory z axis (zL), parallel to B₀. The perpendicular mode is by far the most common.

Beam experiments:
For experiments with a microwave (or THz) beam, use Exp.mwMode = {k, pol}. k specifies the propagation direction in the lab frame, in one of three possible ways: (i) a letter code for the direction, e.g. 'y', 'z', 'xy'; (ii) a 3-element cartesian vector, for example [0;1;0] specifies the lab y axis; (iii) two polar angles [phi_k, theta_k] that specify the orientation. theta_k is the angle between the microwave propagation direction and the lab z axis, and phi_k is the angle between the lab x axis and the projection of the propagation vector onto the lab xy plane. For example, [pi/2, pi/2] gives the lab y axis.

For linearly polarized mw irradiation, additionally provide pol, the polarization angle of the radiation, in radians. To calculate the microwave propagation direction nk and the B1 direction nB1 from k and pol, use

k = 'y';   % propagation along y lab axis
pol = -pi/2; % B1 along x lab axis
[phi,theta] = vec2ang(k);  % convert to angles
[nB1,~,nk]  = erot([phi,theta,pol],'rows')

For unpolarized excitation, set pol='unpolarized'. For circularly polarized radiation, set pol='circular+' or pol='circular-', depending on the sense of rotation.

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

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

SampleFrame is only used for crystals and partially ordered samples (films etc). It is ignored for disordered samples (powders and frozen solutions).

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 (such as 'P212121' or 'Ia-3d'), or the symbol for the point subgroup of the space group (in either Schönflies or Hermann-Mauguin notation, such as 'D2h' or 'mmm').

Exp.CrystalSymmetry = 11;       % space group number (between 1 and 230)
Exp.CrystalSymmetry = 'P21/c';  % space group symbol
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 generated and included in the calculation. If CrystalSymmetry is not given, space group 1 (P1, point group C1, one site per unit cell) is assumed.

Specifies the sample rotation as {nRot_L,rho}. nRot_L fixed in the laboratory frame. The axis is represented in lab-frame coordinates, and it does not need to be normalized. rho is the rotation angle (in radians), or list of rotation angles, around the rotation axis. This rotation is applied to sample starting from the orientation given in Exp.SampleFrame.

rho = pi/4;
Exp.SampleRotation = {[0;1;0],rho};       % rotation around lab y axis (yL)
Exp.SampleRotation = {'z',rho};           % rotation around lab z axis (zL)
Exp.SampleRotation = {'x',rho};           % rotation around lab x axis (xL)

rho = deg2rad(0:30:180);
Exp.SampleRotation = {'x',rho};           % a series of rotations

The three Euler angles, in radians, for the transformation of the sample/crystal frame to the molecular frame. Use this field when specifying a crystal containing spin systems that are tilted with respect to the crystal frame. E.g. Exp.MolFrame=[0,pi/4,0] tilts the x and z axis of the spin system's molecular frame relative to those of the crystal frame while keeping the y axes aligned.

If a number of function handle is given in this field, it specifies the orientational distribution of the paramagnetic molecules in the sample. If not given or set to zero, the distribution is isotropic, i.e. all orientations occur with the same probability.

If a number is given, the orientational distribution is non-isotropic and computed as P(α,β,γ) = exp(-U(β)) with U(β) = -λ(3 cos²β - 1)/2, where λ is the number specified in Exp.Ordering. The angles α, β and γ (in radians) are the Euler angles that transform the sample frame to the molecular frame. In particular, β is the angle between the sample-frame z axis and the molecular-frame z axis.

Typical values for λ are between about -10 and +10. For negative values, the orientational distribution function P(α,β,γ) is maximal at β = 90° (preferential orientation of molecular z axis in the sample xy plane), for positive values it is maximal at β = 0° and β = 180° (preferential alignment of the molecular z axis with the sample z axis). The larger the magnitude of λ, the sharper the distributions.

To plot the orientational distribution for a given value of λ, use

lambda = 2;
beta = linspace(0,pi,1001);
U = -lambda*plegendre(2,0,cos(beta));
%U = -lambda*(3*cos(beta).^2-1)/2;  % equivalent
P = exp(-U);
plot(beta*180/pi,P);

If Exp.Ordering is a function handle, pepper calls the function to obtain the orientational distribution. The function must accept three input arguments (alpha, beta and gamma), in radians. The function must accept arrays and return an array P containing probabilities for each orientation, that is P(k) is the probability of finding the spin centers with orientation (relative to the sample frame) specified by alpha(k), beta(k) and gamma(k). Here is an example with an anonymous function:

Exp.Ordering = @(alpha,beta,gamma) gaussian(beta,0,15*pi/180);

If the function cannot be written as a single line, define a separate function and provide it via its function handle. For example, here is a function oridist defining a von Mises-Fisher distribution

function P = oridist(alpha,beta,gamma)
gamma0 = 0;
beta0 = pi/4;
kappa = 10;
mu = ang2vec(gamma0,beta0);
x = ang2vec(gamma(:),beta(:));
P = kappa/(4*pi)/sinh(kappa)*exp(kappa*mu'*x);
P = reshape(P,size(beta));
end

Include this in the simulation via Exp.Ordering = @oridist.

When providing a function for Exp.Ordering, make sure that the grid symmetry used in the simulation is the same or lower than the symmetry of the distribution. Otherwise, incorrect spectra are obtained. In cases of doubt, set Opt.GridSymmetry='C1'. This always gives correct results.

If the orientational distribution provided via Exp.Ordering is very narrow, increase the number of knots in Opt.GridSize.

The distributions in Exp.Ordering do not have to be normalized. pepper normalizes it internally. If you want to calculate the normalization factor, use numerical integration with integral3:

lambda = 3;
P = @(alpha,beta,gamma) exp(lambda*(3*cos(beta).^2)/2);
integral3(P,0,2*pi,0,pi,0,2*pi)

ans =
   3.2090e+03

Specifies mode of photoexcitation. If photoexcitation is present, photoselection weights will be calculated and included into the spectral line intensities. For this, the transition dipole moment direction must be provided in Sys.tdm (see here). Possible settings for lightBeam are

• Exp.lightBeam = '' - no photoexcitation (default)
• Exp.lightBeam = 'perpendicular' - polarized light incident along the lab y direction (perpendicular to B0) with E-field along the lab x direction (also perpendicular to B0)
• Exp.lightBeam = 'parallel' - polarized light incident along the lab y direction (perpendicular to B0) with E-field along the lab z direction (parallel to B0)
• Exp.lightBeam = 'unpolarized' - unpolarized light incident along the lab y direction (perpendicular to B0) with E-field uniformly distributed in the lab xz plane
• Exp.lightBeam = {k alpha} - light incident along propagation direction k and polarization angle alpha. There are three ways to specify k: (i) a letter code for the direction, e.g. 'y', 'z', 'xy'; (ii) a 3-element vector, e.g. [0;1;0] specifying the lab y axis; (iii) two angles [phi_k theta_k] that specify the orientation. theta_k is the angle between propagation vector and the lab z axis, and phi_k is the angle between the lab x axis and the projection of the propagation vector onto the lab xy plane. For example, [pi/2 pi/2] gives the lab y axis. alpha is the polarization angle, in radians. To represent an unpolarized beam, set alpha=NaN.

Contribution of isotropically scattered light to photoexcitation (see lightBeam), as a value between 0 and 1. 0 means no isotropic contribution (only direct excitation by polarized or unpolarized beam, as given in lightBeam), 1 indicates 100 percent isotropic contribution (essentially beam direction and polarization have no effect).

Determines the level of theory pepper uses to compute the resonance fields (for field sweeps) or frequencies (for frequency sweeps).

• 'matrix' indicates matrix diagonalization. This method is very reliable and accurate and works for spin systems with any number of spins. All interactions, including quadrupole, are included in the computation.
• 'perturb1' indicates first-order perturbation theory, and 'perturb' or 'perturb2' indicates second-order perturbation theory. These methods are limited to spin systems with one electron spin 1/2 (and possibly some nuclei). In addition, nuclear Zeeman and nuclear quadrupole terms are neglected, and only allowed transitions are computed. For multi-nuclear spin system, cross-nuclear effects are neglected as well. The resulting spectrum is reasonably correct only for small hyperfine couplings (e.g. organic radicals).
• 'hybrid' indicates matrix diagonalization for all the electron spins, and perturbation treatment for all nuclei, using effective nuclear sub-Hamiltonians for each electron spin manifold. This method is advantageous for high-spin systems with significant zero-field splitting, but only small hyperfine couplings. If some nuclei have large hyperfine couplings that need to be treated exactly, they can be specified in the field HybridCoreNuclei.

'matrix' is the method of choice for systems with only a few low-spin nuclei (and any number of electron spins). For spin systems with many nuclei and small hyperfine couplings, simulations using perturbation theory are orders of magnitude faster. 'hybrid' is the method of choice for systems with several large electron spins coupled to several nuclei such as in oligometallic clusters.

List of nuclei to include in the matrix diagonalization when using Opt.Method='hybrid'. If not given, it is set to [], and all nuclei are treated perturbationally.

Here is an example:

Sys.Nucs = '63Cu,14N,1H';

Opt.Method = 'hybrid';
Opt.HybridCoreNuclei = [1];   % 63Cu is treated exactly, 14N and 1H perturbationally
Opt.HybridCoreNuclei = [1 2]; % 63Cu and 14N are treated exactly, 1H perturbationally
Opt.HybridCoreNuclei = [];    % all nuclei are treated perturbationally

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. If 'transitions' is specified, transitions subspectra are returned (for powders only). If 'sites' is specified, site subspectra are returned (for crystals only). If 'orientations' is specified, orientation subspectra are returned (for crystals only).

Determines how much information pepper prints to the command window. If Opt.Verbosity=0, pepper is silent. 1 logs relevant information, 2 gives more details.

Determines the number of orientations (knots) in a powder simulation for which spectra are calculated.

• N1 gives the number of orientations between θ=0° and θ=90° for which spectra are explicitly calculated using the physical theory. Common values for N1 are between 10 (10° increments) and 91 (1° increments). The larger the anisotropy of the spectrum and the narrower the line width relative to the anisotropy, the higher N1 must be to yield smooth powder spectra.
• N2 is the refinement factor for the interpolation of the orientational grid. E.g. if N2=4, then between each pair of computed orientations three additional orientations are calculated by spline interpolation. Values higher than 10 are rarely necessary. If N2 is not given, a default value is used.

Opt.GridSize = 91;       % 1° increments, no interpolation
Opt.GridSize = [46 0];   % 2° increments, no interpolation
Opt.GridSize = [31 6];   % 3° increments, 6-fold interpolation (giving 0.5° increments)

Determines the symmetry of the orientational grid used for the powder simulation. 'Dinfh' corresponds to a line from θ=0° to θ=90° (with φ=0°), 'D2h' to one octant, 'C2h' to two octants, 'Ci' to one hemisphere (four octants), and 'C1' to the full sphere. If not given or '', then pepper determines the required grid symmetry automatically from the spin system parameters. With any other setting, pepper will use the specified grid symmetry, even if it is not optimal or incorrect for the spin system. See also hamsymm.

Determines manually the level pairs which are used in the spectrum calculation. If given, pepper uses them and skips its automatic transition selection scheme. Level pairs are specified in Transitions(k,:) by the level numbers which start with 1 for the lowest-energy level. If 'all' is given, then all transitions are included.

Opt.Transitions = [1 6; 2 5];       % include only transitions 1->6 and 2->5

In crystal simulations, this gives a list of crystal sites to include in the simulation.

If Opt.Sites is empty or not given, all sites are included. If given, it must be a list of site numbers. The number of sites depends on the space group given in Exp.CrystalSymmetry. E.g. the following set limits the simulation to sites 1 and 3 of the 4 magnetically distinct sites in crystal of space group no. 47.

Exp.CrystalSymmetry = 47;  % space group Pmmm
Opt.Sites = [1 3];

In powder simulations, Opt.Sites is ignored.

Specifies the threshold for transition pre-selection. Only transitions with an estimated relative average intensity larger than this number are included. The relative average intensity of the strongest transition is set to 1. The default value for the threshold is 1e-4. The pre-selection is an approximate procedure, and it might miss transitions for complicated spin systems. In these cases, setting it to zero will include all transitions in the simulation. After the intensities of all included transitions are computed, the transition are screened again against this threshold in a post-selection step. If transitions are specified manually in Opt.Transitions, Opt.Threshold is ignored.

With 'on', transition rates, i.e. line intensities, are computed correctly. Allowed transitions will be more intense then quasi-forbidden ones. 'off' simply sets all transition rates of all transitions to 1. Allowed and forbidden transitions will have the same intensity. Be very careful when switching this option to 'off'! The resulting spectra are not correct.

Determines whether the frequency-to-field conversion factor is included in the line intensities of field-swept spectra. 1 indicates yes, 0 indicates no. The factor is the generalized 1/g Aasa-Vänngård factor. This setting is ignored for frequency-swept spectra.

For isotope mixtures, determines which isotopologues to include in the simulation. Any isotopologue with relative abundance smaller than IsoCutoff is excluded. The default value is 1e-4.

The amount of random noise to add to non-zero Hamiltonian matrix elements to break degeneracies. This is needed in some cases, since EasySpin cannot handle systems with degeneracies in a numerically stable fashion. The default value is 1e-10. This means each non-zero Hamiltonian matrix element is multiplied by a random number between (1-1e-10) and (1+1e-10). This is only used when matrix diagonalization is used.

Sys = struct('S',1/2,'g',[1.9 2 2.3]);
Exp = struct('CenterSweep',[325 80],'mwFreq',9.5,'Harmonic',0);
Opt = struct('Verbosity',1);

[x,y1] = pepper(Sys,Exp,Opt);
Sys.lw = 2;
y2 = pepper(Sys,Exp,Opt);
Sys.lw = 0;
Sys.HStrain = [170 40 50];
y3 = pepper(Sys,Exp,Opt);

plot(x,y1/sum(y1),x,y2/sum(y2),x,y3/sum(y3));
legend('no broadening','convolution broadening','H strain');