pepper
Synopsis

Calculation of single-crystal and powder cw EPR spectra (solid-state).

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

See also the user guide on how to use pepper.

Description

pepper calculates cw EPR spectra for powders, frozen solutions and single crystals.

There are up to three possible output arguments

There are three arguments to the function, the last one optional. They are similar to those of the function resfields.

Sys is a spin system structure containing spin Hamiltonian parameters and line broadening parameters.

For simulating a multi-component mixture, Sys should be a cell array of spin systems, e.g. {Sys1,Sys2} for a two-component mixture. Each of the component spin systems should have a field weight that specifies the weight of the corresponding component in the final spectrum.

Exp contains experimental parameters such as the microwave frequency, the magnetic field range and temperature. Here is a full list of its fields:

mwFreq Spectrometer frequency, in GHz. E.g. Exp.mwFreq = 9.5; for X band.
CenterSweep 2-element vector [center sweep]
Contains center field and sweep width of the external magnetic field sweep range. Values should be in mT, e.g. Exp.CenterSweep=[310 100].

The magnetic field sweep range can be specified either in CenterSweep or in Range. If both are given, CenterSweep has precedence.

Range 2-element vector [loField hiField]
Contains lower and upper limit of the external magnetic field sweep range. Values should be in mT, e.g. Exp.Range=[260 360].

The magnetic field sweep range can be specified either in CenterSweep or in Range. If both are given, CenterSweep has precedence.

nPoints Number of points on the magnetic field abscissa axis. If not given, the default is 1024.
Temperature scalar (default inf) or vector

This field specifies populations for the states of the spin system, either directly or via a temperature.

Thermal equilibium:
Temperature is the temperature of the spin system in the EPR experiment, in Kelvin. If given, Boltzmann populations are computed and included in the EPR line intensities. E.g., Temperature = 298 corresponds to room temperature. If not given (or set to inf), all transitions are assumed to have equal polarizations.

Non-equilibrium populations:
Temperature can also be used to specify non-equilibrium populations. For a spin system with N electron states (e.g. 4 for a biradical), it must be a vector with N elements giving the populations of the zero-field electron states, from lowest to highest in energy.
E.g., if Temperature = [0.85 0.95 1.2] for an S=1 system, the population of the lowest-energy zero-field state is 0.85, and the highest-energy zero-field state has a population of 1.2. The populations of all nuclear sublevels within an electron spin manifold are assumed to be equal.

Harmonic 0, 1 (default) or 2
Harmonic of the detection. 1 and 2 specify the first and the second derivative of the absorption spectrum, respectively. 0 returns the absorption spectrum without differentiation. To explicitly include the effect of field modulation, use Exp.ModAmp.
Mode 'perpendicular' (default) or 'parallel'
Relative orientation of the microwave field B1 with respect to the static magnetic field B0. 'perpendicular' means [eqn]. 'parallel' means that [eqn].
ModAmp Modulation amplitude (peak-to-peak), in mT.
mwPhase The reference microwave phase, in radians. 0 is pure absorption (default value), and pi/2 is pure dispersion. mwPhase is used only if the convolutional broadening given in Sys.lw or Sys.lwpp has a Lorentzian component.

mwPhase allows you to include absorption/dispersion admixture in the simulation.

Orientations 3xN array of real
Specifies single-crystal orientations for which the EPR spectrum should be computed. Each column of Orientation contains the three Euler rotation angles [phi;theta;chi] of the magnetic field in a molecule fixed frame. If Orientation is empty or not specified, the full powder spectrum is computed.
Exp.Orientations = [0;0;0];              % crystal with z axis aligned with B0
Exp.Orientations = [0;pi/2;0];           % crystal with z axis perpendicular to B0
Exp.Orientations = [0 0 0;0 pi/2 0].';   % two orientations
Exp.Orientations = [];                   % powder
CrystalSymmetry

Specifies the symmetry of the crystal, if single-crystal spectra to be simulated (that is, if Exp.Orientations is specified). The crystal symmetry can be either the number of the space group (between 1 and 230), the symbol of the space group (e.g. 'P21212' or the symbol for the point group (e.g. 'C2h' or '2/m').

Exp.CrystalSymmetry = 'P21/m'; % 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, pepper automatically computes the spectra of all symmetry-related sites in the crystal. If CrystalSymmetry is not given, pepper assumes space group 1 (P1, point group C1), which has only one site per unit cell.

Ordering scalar (default: zero) or function handle

If a number 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 it is given, the orientational distribution is non-isotropic and computed according to the formula P(θ) = exp(λ(3 cos2θ - 1)/2), where θ is the angle between the molecular z axis and the static magnetic field, and λ is the number specified in Exp.Ordering.

Typical values for λ are between about -20 and +20. For negative values, the orientational distribution function is maximum at θ = 90°, for positive values at θ = 0° and θ = 180°. The larger the magnitude of λ, the sharper the distributions.

To plot a distribution depending on λ, use

lambda = 5;
theta = linspace(0,pi,1001);
p = exp(lambda*plegendre(2,0,cos(theta)));
plot(theta*180/pi,p);

If Exp.Ordering is a function handle, pepper will use the user-supplied function to obtain the orientational distribution. It calls the function with two vector arguments, phi and theta (in radians). The function must return a vector P containing probabilities for each orientation, that is P(k) is the probability of finding the paramagnetic molecules with orientation specified by phi(k) and theta(k). Here is an example with an anonymous function:

Exp.Ordering = @(phi,theta) gaussian(theta,0,15/180*pi);

Of course, the function can also be written and stored in a separate file, e.g. myori.m. Then use Exp.Ordering = @myori.

When using a custom orientational distribution, make sure that the symmetry used in the simulation corresponds to the symmetry of the distribution. If the distribution is very narrow, increase the number of knots in the options structure.

mwFreq and Range have to be provided by the user, all other fields are optional and have default values.

The structure Opt collects computational parameters. Opt need not be specified, in which case default values for all fields are used. The field names and their possible values are listed below.

Verbosity Determines how much information pepper prints to the screen. If Opt.Verbosity=0, pepper is completely silent. 1 logs relevant information, 2 gives more details.
Output 'summed' (default) or 'separate'
Determines in what form the spectrum is returned. If set to 'separate', pepper returns one spectrum for each transition in a matrix spec. The transition spectra are along the rows. spec(k,:) is the spectrum of transition k. If 'summed' is specified, the total spectrum is returned in spec as a vector.
nKnots [N1] or [N1 N2]
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 explicitely 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 linewidth, 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.nKnots = 91;       % 1° increments, no interpolation
Opt.nKnots = [46 0];   % 2° increments, no interpolation
Opt.nKnots = [31 6];   % 3° increments, 6-fold interpolation (giving 0.5° increments)
Symmetry 'auto' (default), 'Dinfh', 'D2h', 'C2h' or 'Ci'
Determines the symmetry used for the powder simulation. Based on this the set of orientations for which spectra are computed is chosen. 'Dinfh' corresponds to a line from θ=0° to &theta=90° (with φ=0°), 'D2h' to one octant, 'C2h' to two octants, and 'Ci' to one hemisphere (four octants). auto is the default, meaning that pepper determines the correct symmetry automatically from the given spin system. With any other setting, pepper is forced into using the specified symmetry, even if it is incorrect for the spin system. See also symm.
Transitions mx2 vector of integers
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.
Threshold Specifies the threshold for pepper's transition pre-selection. Any transition with an estimated relative average amplitude less than this number is not included in the calculation. The relative average amplitude of the strongest transition is 1, the default 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 computation.
Intensity 'on' (default) or 'off'
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.
Method 'matrix' (default), 'perturb', 'perturb1', 'perturb2', 'hybrid'

Determines how pepper computes the resonance fields.

  • '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 ares 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. 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.
'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.

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

Spectra are calculated over a triangular orientational grid using resfields to obtain the resonance field positions and amplitudes. For each orientation line positions, and possibly widths and intensities, are evaluated.

This gridded data is then interpolated with cubic splines in a combined 1D/2D approach. Resampling of the spline surface gives much quicker many more position/intensity/width data than quantum-mechanical calculation.

Finally, the refined data are projected onto the magnetic field axis using a Delaunay triangulation of the resampled spline surfaces. Linear interpolative projection of these triangles yields an extremely smooth spectrum with very low powder simulation noise. In the case of full anisotropic width treatment, a simple sum-up of Gaussian line shapes is used instead of the projection.

Apart from the main steps above, there is an automatic transition selection, which works along the same line as the overall algorithm, except that its results are only used for determining which level pairs possibly contribute to the spectrum.

For line width calculations, Gaussian distributions are assumed both in the magnetic field and the frequency dimension. The overall line width for a given orientation is

[eqn]

where [eqn] is the residual line width specified in Sys.HStrain, [eqn] is the line width due to correlated g-A strain (Sys.gStrain and Sys.AStrain), and [eqn] the width arising from D-E strain (Sys.DStrain).

Although quite robust and general, pepper still has some limitations.

Examples

As an illustration, we explore the influence of various pepper options on the zeroth-harmonic (DC) spectrum of a simple orthorhombic system. First the spin system, the experiment at X-band and some options are defined. An anisotropic line width is included in the spin system.

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);

Next we compute spectra for some combinations of broadening parameters.

[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);

The final plot reveals the differences between the spectra.

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

References which contain concepts, formulas and algorithms directly used in the function are listed below.

See also

eigfields, esfit, garlic, resfields, salt, sham