This tutorial explains how to simulate solid-state cw EPR spectra of powders or single crystals using EasySpin. It is assumed that you are familiar with the basic syntax of MATLAB, esp. with structures.
It contains the following topics:pepperSolid-state cw EPR spectra of powders and single crystals are computed by the EasySpin function pepper. It can be called with two or three parameters and can return both field axis and spectrum.
pepper(Sys,Exp); % plots the spectrum [Field,Spec] = pepper(Sys,Exp); % returns the field axis and the spectrum [Field,Spec] = pepper(Sys,Exp,Opt); % additional simulation options in Opt
Don't forget the ; (semicolon) at the end of the line to suppress output to the screen.
The first argument Sys tells pepper all about the
spin system, and the second argument Exp gives the experimental
parameters. The third, optional, argument Opt contains settings
concerning the simulation itself, like the number of orientations for powder simulations.
The outputs Field and Spec are arrays containing the magnetic field values
and the spectrum, respectively. If no output is requested, pepper simple plots the spectrum.
If the outputs are requested, pepper does not plot the spectrum, but you can plot it
yourself using
plot(Field,Spec);
Setting up a simulation and running it takes only a few lines of code. A very simple one would be
Sys.g = [2 2.1]; Sys.lwpp = 0.5; Exp.mwFreq = 9.5; pepper(Sys,Exp);
This simulates and plots the 9.5 GHz EPR spectrum of an S=1/2 system with an axial g tensor. Copy and paste the code above to your MATLAB command window to see the graph.
The first input argument to pepper is a structure specifying
the spin system.
It contains fields for the electron spin(s), the nuclear spins,
and the various interaction matrices like g and hyperfine tensor.
pepper automatically assumes S=1/2 for the spin quantum number.
For systems with higher spin or more than one unpaired electron, the spin quantum number
should be given in the field S.
Sys.S = 1; % a triplet state Sys.S = 5/2; % for e.g. high-spin Mn2+ or high-spin Fe3+ Sys.S = [1/2, 1/2]; % for a biradical
The field g contains the principal values of the g tensor(s). A simple orthorhombic
S=1/2 system (e.g., a low-spin Fe3+) is
Sys.g = [1.8, 2, 2.1];
Nuclear spins are included by specifying Nucs (comma-separated
list of nuclei) and A (array of hyperfine tensor principal values,
in MHz).
Sys.Nucs = '2H'; % one 2H (deuterium) nucleus Sys.A = [-1,-1,2]*4.2; % hyperfine principal values in MHz
If the A tensor is tilted with respect to the molecular frame, the
tilt angles can be provided via the
field Apa (standing for "A principal
angles")
Sys.Apa = [0 30 0]*pi/180; % [alpha beta gamma] in radians
The zero-field splitting is specified in the D field, in units of
MHz. There are several different input possibilities:
Sys.D = 120; % D = 120 MHz, E = 0 Sys.D = [120 15]; % D = 120 MHz, E = 15 MHz Sys.D = [-25,-55,80]; % principal values of D tensor, in MHz
D and E are related to the principal values of the D tensor (see reference page on the zero-field splitting).
Details about all the spin Hamilton parameters can be found on the spin Hamiltonian reference page. It is also possible to include several electron spins. Refer to the page about spin system structures for details.
No cw EPR spectrum is infinitely sharp. Lines are usually broadened due
to several reasons. pepper provides means to include several
line broadening models in a simulation.
The simplest way to include line broadening is to convolute a stick spectrum
with a (Gaussian or Lorentzian) lineshape after the end of the simulation.
Such a convolution broadening is specified in the spin system field lwpp.
Sys.lwpp = 0.5; % Gaussian broadening of 0.5 mT PP Sys.lwpp = [0 2]; % Lorentzian broadening of 2 mT PP Sys.lwpp = [1 2]; % Gaussian PP of 1mT + Lorentzian PP of 2 mT
The line width is in mT and refers to peak-to-peak (PP) line widhts.
Instead, FWHM (full width at half height) line widths can be provided in the field
Sys.lw.
Sys.lw = 0.5; % Gaussian broadening of 0.5 mT FWHM Sys.lw = [0 2]; % Lorentzian broadening of 2 mT FWHM
For details about line shapes and conversion formulas to/from FWHM and peak-to-peak widths, see the page on line shapes.
Physically, there are several origins for line broadening. Large contribution to broadening often comes from unresolved hyperfine couplings and from distributions in the various magnetic parameters lige g, A and D that result from structural variations from one paramagnetic center to the next.
To include effects from unresolved hyperfine couplings, an orientation-dependent
phenomenological broadening can be specified in HStrain:
Sys.HStrain = [50 50 87]; % [along x, along y, along z], in MHz
Distributions in magnetic parameters are called strains. g and A strains are given in similar fields:
Sys.gStrain = [0.01 0.02 0.005]; Sys.AStrain = [10 10 30]; % in MHz
The three values in gStrain are the FWHM parameters of the Gaussian
distributions of the respective g principal values given in Sys.g.
AStrain is the same for the A tensor. The g and A strains are
correlated.
Distributions of the D tensor values can be given in DStrain, where
the first value is the width of the (scalar) D distribution, and the second is
the width for the E distribution.
All these broadening parameters can be combined. However, usually a modelling
of the broadening with lwpp or HStrain is absolutely
sufficient.
All experimental settings are given in the second input argument Exp.
Just as the spin system, Exp is a structure containing several
fields.
For every simulation, the spectrometer frequency has to
be given in the field mwFreq in units of GHz (Gigahertz).
Exp.mwFreq = 9.754; % in GHz
There are two ways to specify the magnetic field sweep range.
Exp.CenterSweep = [340 80]; % in mT Exp.Range = [300 380]; % in mT
Either the
center and the sweep width (in mT) are given in Exp.CenterSweep,
or the lower and upper limit of the sweep range (again in mT) are given in
Exp.Range.
In many cw EPR spectrometers, the field range is
specified using a center field and a sweep width, so Exp.CenterSweep
is the more natural choice.
Exp.CenterSweep and Exp.Range are only optional.
If both are omitted, pepper tries to automatically determine a
field range large enough to accomodate the full spectrum.
If, on the other hand, both Exp.CenterSweep and Exp.Range are given, pepper takes
the values given in Exp.CenterSweep and ignores those
in Exp.Range.
By default, pepper computes a 1024-point spectrum. However, the number of
points can be changed manually to a different value, e.g.,
Exp.nPoints = 5001;
By default, pepper computes the first-harmonic absorption spectrum, i.e. the first
derivative of the absorption spectrum. By changing Exp.Harmonic, the absorption spectrum
directly or the second-harmonic (second derivative) of it can be requested.
Exp.Harmonic = 0; % absorption spectrum, direct detection Exp.Harmonic = 1; % first harmonic (default) Exp.Harmonic = 2; % second harmonic
If you want to include effects of the modulation amplitude and the time constant, apply the functions fieldmod and rcfilt to the simulated spectrum.
For more advanced spectral simulations, pepper offers more
configuration possibilities in the experimental parameter structure Exp.
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. pepper can simulate both
types of spectra:
Exp.Detection = 'perpendicular'; % perpendicular detection (default) Exp.Detection = 'parallel'; % parallel detection
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, pepper will include the relevant polarization
factors resulting from a thermal equilibrium population of the energy levels.
If not specified otherwise, pepper computes a powder spectrum. But
it is as well straightforward to simulate spectra for a single crystal. The
orientation (or orientations if more than one) of the single crystal can be provided in the
experiment structure field Exp.Orientations. This field should
contain the tilt angles between molecular and laboratory frame
(right-handed coordinate system with z along the static
field and x along the microwave magnetic field), one set of three angles per
column.
For a crystal with its molecular frame aligned with the laboratory frame, the setting is
Exp.Orientations = [0;0;0];
If you need more than one crystal at the same time, then just specify more than one orientation.
Exp.Orientations = [0 0 0;0 pi/4 0].';
In many crystals, there are several symmetry-related sites with identical paramagnetic centers
differing only in their orientation in the crystal. You can tell pepper
about this by providing the crystal symmetry in the field Exp.CrystalSymmetry, e.g.
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
With the crystal symmetry given, pepper not only computes the spectrum for the orientation
given in Exp.Orientations, but also for all symmetry-related sites.
If Exp.Orientations set to [] (an empty array),
pepper simulates the powder spectrum.
The third input argument to pepper contains simulation options.
All of them have reasonable default values, but sometimes it might be necessary
to change one of them. In the following the most important ones are presented.
If you want pepper to print information about the simulation to
the command window during the computation, use
Options.Verbosity = 1;
'Verbosity' tells pepper how much of progress
information to show in the command window. 0 (the default)
suppresses all output, 1 is normal intormation, and 2 prints more information,
relevant only for debugging.
Another useful option is nKnots, which determines how many
orientations pepper will include in the simulation of a powder
spectrum. If this value is too low, the spectrum shape contains ripples.
nKnots is the number of orientations between the z axis and
the x axis (between theta = 0 and theta = 90 degrees).
Options.nKnots = 31; % corresponds to 3-degree increments
The higher nKnots, the finer the orientational grid.
The default value of 19 (5-degree increments) is appropriate for most systems.
A value larger than 181 (0.5-degree increments) is rarely needed.
After having computed the spectrum for a number of orientations specified
by nKnots, the simulation function interpolates these spectra for
additional orientations before summing up all spectra. This interpolative
refinement can be configured with a second number in nKnots.
nKnots = [19 4] means that pepper interpolates
additional 4 spectra between two adjacent orientations evaluated.
Options.nKnots = [19 10]; % massive interpolation Options.nKnots = [19 0]; % no interpolation
The option Output can be used to determine the form in which
pepper returns the spectral data.
% single crystal: orientations separately % powders: transitions separately Options.Output = 'separate'; % sum over all orientations and transitions Options.Output = 'summed';
There are more option fields available. For details, see the documentation page on pepper.
Simulations of spectra from systems with more than one nucleus can be
very time-consuming. To accelerate such computations, pepper provides
two alternatives: (1) nuclei with small hyperfine couplings can be treated by
first-order perturbation theory, still treating electrons and nuclei with large
hyperfine couplings with matrix diagonalization; or (2) the entire system can be treated by first- or second-order
perturbation theory. Let's look at the two possibilities in turn.
The two relevant simulation options are Opt.Method ('matrix' or 'perturb1')
and Opt.Perturb (0 or 1).
As an example, we look at the simulation of the spectrum of Cu2+ porphyrin
Sys.S = 1/2; Sys.g = [2 2.2]; Sys = nucspinadd(Sys,'63Cu',[50 500]); A = [20 30]; Sys = nucspinadd(Sys,'14N',A); Sys = nucspinadd(Sys,'14N',A); Sys = nucspinadd(Sys,'14N',A); Sys = nucspinadd(Sys,'14N',A); Sys.lwpp = 0.4;
With matrix diagonalization (Opt.Method='matrix', which is the default), the simulation needs several hours. With the perturbational
methods, the simulation is orders of magnitude faster, but less accurate.
First, let's run the simulation with selective perturbational treatment of nuclei
with small hyperfine couplings. To enable this perturbational treatment of super-hyperfine
(SHF) nuclei, set the field Perturb in the options structure to 1.
Exp.mwFreq = 9.5; Exp.Range = [270 360]; Opt.nKnots = 31; Opt.Perturb = 1; [x,y] = pepper(Sys,Exp,Opt);
pepper will automatically determine those nuclei where first-order
perturbation theory can be applied without introducing noticeable errors in the
simulated spectrum. The simulation takes a few seconds only.
The second perturbational method, which uses second-order perturbation theory for the entire spin system, gives a simulation that is much faster, but also much less accurate.
Opt.Method = 'perturb2'; pepper(Sys,Exp,Opt);
This tells pepper to use second-order perturbation theory. The simulation
needs only a small fraction of a second. If you want to compare the full matrix diagonalization
to the perturbation simulation, use
Opt.Method = 'matrix'; [x,y1] = pepper(Sys,Exp,Opt); Opt.Method = 'perturb2'; [x,y2] = pepper(Sys,Exp,Opt); plot(x,y1/max(y1),'r',x,y2/max(y2),'b');
pepper can handle both thermal equilibrium and non-equilibrium
populations. Both are specified in the field Temperature of the
experimental settings structure.
For thermal equilibrium, just give the temperature in Kelvin:
Exp.Temperature = 77; % 77K, boiling point of liquid nitrogen
For non-equilibirum populations, Temperature must be a
vector. If the spin systems contains N electron states, then
this vector must contain N elements, each specifying the population
of one of the electron states at zero field, sorted according to
their energy from lowest to highest.
E.g., an organic triplet with S=1 and I=1 has 3 electron states, each further split into three sublevels by the coupling to the nuclear spin. The population vector in this case should contain three elements:
Exp.Temperature = [0.6 0.8 1.1]; % highest state is most populated
This specifies that all the sublevels of the lowest zero-field electron states
have a population of 0.6, etc. The sublevels of the highest-energy zero-field
electron state have a population of 1.1. The populations don't have to be normalized,
pepper takes care about that.
To compute the state populations for a non-zero field state,
pepper decomposes it into a linear combination of zero-field
states and combines the zero-field populations using the resulting linear
combination coefficients.
A simple example of a non-equilibrium triplet system is
Sys.S = 1; Sys.g = 2; Sys.lw = 0.2; Sys.D = 100; Exp.mwFreq = 9.5; Exp.Range = [320 360]; Exp.Harmonic = 0; Exp.Temperature = [0.5 0.6 0.9]; pepper(Sys,Exp);
In powders and frozen solutions (disordered systems), paramagnetic molecules are randomly oriented in the sense that each orientation occurs with equal probability. In other situations, like in polymers, biomembranes or liquid crystals, the paramagnetic molecules may be partially aligned or ordered, so that some orientations are more probable than others. As a result, the spectra of such partially ordered systems are different from those of powders and frozen solutions.
pepper can include partial ordering in the spectral simulation. For this, set a value
(different from zero) in the experiment structure field Exp.Ordering. It is a number which specifies the nature and
the degree of the ordering.
If it is zero, no ordering is used. Positive values mean that the molecules are partially aligned so that the magnetic field vector is mostly approximately parallel to the molecular z axis. Negative values specify preferential orientation away from the z axis, that is, orientations with the magnetic field vector in the xy plane are more probable.
Exp.Ordering = 0; % all orientations equally populated Exp.Ordering = -1; % slightly preferential orientation in the xy plane Exp.Ordering = 10; % strongly aligned along the z axis
The ordering function used is very simple (see the pepper documentation for more information), but sufficient for most cases. Here is a simulation of a sample where the molecules are preferentially oriented so that the molecular z axis is close to the magnetic field vector:
Sys.g = [2 2 2.2]; Sys.lwpp = 1; Exp.mwFreq = 9.5; Exp.Ordering = +2; pepper(Sys,Exp);
You can also define your own custom orientational distribution in a separate function and supply to
pepper as a function handle in Exp.Ordering. See the pepper documentation for details. Written as an anonymous function, the built-in orientational distribution is equivalent to
Exp.Ordering = @(phi,theta) exp(lambda*plegendre(2,0,cos(theta)));
where the EasySpin function plegendre returns the associated Legendre polynomial and lambda corresponds to the number given in Exp.Ordering.