Simulating solid-state cw EPR spectra

This user guide explains how to simulate solid-state cw EPR spectra of powders, frozen solutions, partially ordered samples and single crystals using EasySpin's function pepper. It is assumed that you are familiar with the basics of MATLAB, esp. with structures.

It contains the following topics:

There are the following advanced topics:

Running the simulation

Solid-state cw EPR spectra 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
[B,spec] = pepper(Sys,Exp);       % returns the field axis and the spectrum
[B,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 (sample properties and spectrometer settings). The third, optional, argument Opt contains settings concerning the simulation itself, like the number of orientations for powder simulations.

The outputs B and spec are arrays containing the magnetic field values and the spectrum, respectively. If no output is requested, pepper plots the spectrum. If the outputs are requested, pepper does not plot the spectrum, but you can plot it yourself using 

plot(B,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;    % mT
Exp.mwFreq = 9.5;  % GHz
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 plot.

The spin system

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 parameters like g and hyperfine tensors.

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(s) should be given in the field S. 

Sys.S = 1;           % a triplet
Sys.S = 5/2;         % a sexted, e.g. high-spin Mn(II) or high-spin Fe(III)
Sys.S = [1/2, 1/2];  % a biradical

The field g contains the principal values of the g tensor(s). A simple rhombic 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 AFrame 

Sys.AFrame = [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 MHz
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.

Broadenings

No cw EPR spectrum is infinitely sharp. Lines are usually broadened for 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) line shape after the end of the simulation. Such a convolutional 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 1 mT + Lorentzian PP of 2 mT

The line width is in mT and refers to peak-to-peak (PP) line widths in the first-harmonic EPR spectrum. FWHM (full width at half height) line widths can be provided instead 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 between FWHM and PP widths, see the page on line shapes.

Physically, there are several possible origins for line broadening. Large contributions to broadening often comes from unresolved hyperfine couplings and from distributions in the various magnetic parameters like g, A and D that result from structural variations from one paramagnetic center to the next in the sample.

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 for uncorrelated 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 can be set to be correlated using Sys.gAStainCorr.

Distributions of the D tensor values can be given in DStrain, where the first value is the FWHM of a Gaussian distribution of D, and the second is the width for a Gaussian distribution of E.

All these broadening parameters can be combined. However, usually a modelling of the broadening with lwpp or HStrain is sufficient.

Basic experimental settings

The second input argument, Exp, collects experimental settings. This includes additional information about the sample as well as spectrometer settings. Just as the spin system, Exp is a structure containing several fields.

Microwave frequency. To simulate an EPR spectrum, Easyspin needs at a minimum the spectrometer frequency. Put it into Exp.mwFreq, in units of GHz. 

Exp.mwFreq = 9.385;  % X-band
Exp.mwFreq = 34.9;   % Q-band

Field range. There are two ways to provide the magnetic field sweep range. Either give the center field and the sweep width (in mT) in Exp.CenterSweep, or specify the lower and upper limit of the sweep range (again in mT) in Exp.Range. 

Exp.CenterSweep = [340 80]; % in mT
Exp.Range = [300 380];      % in mT

If Exp.CenterSweep and Exp.Range are omitted, EasySpin tries to determine a field range large enough to accommodate the full EPR spectrum. This automatic ranging works for many simple systems, but fails in more complicated situations. EasySpin will issue an error when it fails.

Points. By default, pepper computes a 1024-point spectrum. However, you can change the number of points to a different value using 

Exp.nPoints = 5001;

Harmonic. By default, EasySpin computes the first-harmonic absorption spectrum, i.e. the first derivative of the absorption spectrum. By changing Exp.Harmonic, you can request the absorption spectrum directly or the second-harmonic (second derivative) of it. 

Exp.Harmonic = 0;  % absorption spectrum, direct detection
Exp.Harmonic = 1;  % first harmonic (default)
Exp.Harmonic = 2;  % second harmonic

Modulation amplitude. If you want to include effects of field modulation like overmodulation, use Exp.ModAmp. 

Exp.ModAmp = 0.2;   % 0.2 mT (2 G) modulation amplitude, peak-to-peak

Time constant. To include the effect of the time constant of an RC filter, apply the function rcfilt to the simulated spectrum.

More experimental settings

For more advanced spectral simulations, EasySpin offers more possibilities in the experimental parameter structure Exp.

Mode. 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 magnetic field. Some resonators can operate in parallel mode, where the microwave field is parallel to the static one. EasySpin can simulate both types of spectra: 

Exp.mwMode = 'perpendicular'; % perpendicular mode (default)
Exp.mwMode = 'parallel';      % parallel mode

Temperature. 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, EasySpin includes the relevant polarization factors resulting from a thermal equilibrium population of the energy levels. For S=1/2 systems, including the temperature typically does not affect the spectral shape. However, temperature strongly affects the spectral shape for high-spin systems with large zero-field splittings, and for coupled spin systems with significant exchange couplings.

Microwave phase. Occasionally, the EPR absorption signal has a small admixture of the dispersion signal. This happens for example when the microwave phase in the reference arm is not absolutely correctly adjusted. EasySpin can mix dispersion with absorption if a Lorentzian broadening is given: 

Sys.lwpp = [0.2 0.01];           % Lorentzian broadening (2nd number) required

Exp.mwPhase = 0;                 % pure absorption
Exp.mwPhase = pi/2;              % pure dispersion
Exp.mwPhase = 3*pi/180;          % 3 degrees dispersion admixed to absorption
Powders and crystals

If not specified otherwise, pepper computes a powder spectrum, i.e. it assumes that the orientations of the spin centers are randomly but uniformly distributed in the sample. To simulate the spectrum for a single crystal, use the experiment structure field Exp.MolFrame to provide the orientation of the spin center within the crystal, Exp.CrystalSymmetry to provide the symmetry of the crystal, and Exp.SampleFrame to specify how the crystal is oriented in the spectrometer. For details see the dedicated page on crystal simulations.

Exp.MolFrame contains the tilt angles that transform the sample/crystal frame to the molecular frame. 

Exp.MolFrame = [0 0 0];          % molecular frame aligned with sample frame
Exp.MolFrame = [0 30 0]*pi/180;  % molecular frame titled

In most crystals, there is a set of symmetry-related sites with identical paramagnetic centers differing only in their orientation in the crystal. This depends on the space group of the crystal. Provide information about this 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, EasySpin includes all symmetry-related sites into the simulated spectrum.

Finally, use Exp.SampleFrame to specify how the crystal is oriented relative to the laboratory frame. For a crystal sample with its frame aligned with the laboratory frame, the setting is 

Exp.SampleFrame = [0 0 0];   % sample/crystal frame aligned with lab frame

If Exp.SampleFrame is missing or set to [] (an empty array), pepper simulates the powder spectrum.

Simulation options

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 

Opt.Verbosity = 1;

Opt.Verbosity tells pepper how much of progress information to show in the command window. 0 (the default) suppresses all output, 1 is normal information, and 2 prints more information, relevant only for debugging.

Another useful option is Opt.GridSize, which determines how many orientations pepper includes in the simulation of a powder spectrum. If this value is too low, the spectrum shape contains ripples. Opt.GridSize is the number of orientations between the z axis and the x axis (between theta = 0 and theta = 90 degrees). 

Opt.GridSize = 31;   % corresponds to 3-degree increments

The higher Opt.GridSize, the finer the orientational grid. The default value of 19 (5-degree increments) is appropriate for many 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 Opt.GridSize, 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 Opt.GridSize. Opt.GridSize = [19 4] means that pepper interpolates additional 4 spectra between two adjacent orientations evaluated. 

Opt.GridSize = [19 10];  % massive interpolation
Opt.GridSize = [19 0];   % no interpolation

The option Opt.separate can be used to obtain a set of individual transition spectra instead of the total spectrum: 

Opt.separate = 'transitions';  % separate subspectra for all transitions
Opt.separate = '';             % one overall spectrum

There are more option fields available. For details, see the documentation page on pepper.

Multiple components

Often, an EPR spectrum shows a mixture of spin species. To simulate these spectra, each of the component spectra has to be simulated and added with the appropriate weight (depending on spin concentration) to the total spectrum.

This can be done automatically by pepper. Just provide the component spin systems with their weights as a cell array (in curly braces) to pepper. For example, here is the simulation of a very simple two-component mixture with 2:1 ratio of spin concentrations. 

Sys1.g = 2;
Sys1.lwpp = 1;   % mT
Sys2.g = 2.1;
Sys2.lwpp = 0.8;  % mT

Sys1.weight = 2;
Sys2.weight = 1;

Exp.mwFreq = 9.5;  % GHz
Exp.Range = [300 360];  % mT

pepper({Sys1,Sys2},Exp);

You don't have to specify Sys.weight - if it's absent it is assumed to be 1. These weights are absolute, i.e. a simulation with Sys.weight=20 yields a spectrum that is 10 times more intense than the one obtained with Sys.weight=2. There is no limit to the number of components in a simulation.

Systems with several nuclei

pepper uses matrix diagonalization as the default method for simulating spectra. For systems with several nuclei this can be very time-consuming. To accelerate such computations, pepper provides first- and second-order perturbation theory as an alternative methods. The relevant simulation option that tells EasySpin about is Opt.Method.

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];   % MHz
Sys = nucspinadd(Sys,'14N',A);
Sys = nucspinadd(Sys,'14N',A);
Sys = nucspinadd(Sys,'14N',A);
Sys = nucspinadd(Sys,'14N',A);
Sys.lwpp = 0.5;  % mT

With matrix diagonalization (Opt.Method='matrix', which is the default), the simulation needs several hours. With second-order perturbation theory (Opt.Method='perturb2'), the simulation is orders of magnitude faster, but potentially less accurate. We can compare the full matrix diagonalization to the perturbation simulation. 

Exp.mwFreq = 9.5;   % GHz
Exp.Range = [260 380];  % mT
Opt.Method = 'matrix';
[x,y1] = pepper(Sys,Exp,Opt);
Opt.Method = 'perturb2';
[x,y2] = pepper(Sys,Exp,Opt);
plot(x,y1,'r',x,y2,'b');
Partially ordered systems

In powders and frozen solutions (disordered systems), paramagnetic molecules are randomly oriented in the sense that any orientation can occur 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. Information about the ordering is provided via Exp.Ordering. There are two ways to specify ordering, either using a single number or by providing an anonymous function. If Exp.Orderingis zero or empty, no ordering is used.

If Exp.Ordering is set to a number, a simple axial built-in ordering function is used. Positive values mean that the molecules are partially aligned so that the molecular z axis is predominantly along the sample z axis. Negative values specify partial ordering such that molecular orientations with the molecular z axis perpendicular to the sample z axis are dominant. 

Exp.Ordering = 0;   % all orientations equally populated
Exp.Ordering = -1;  % slightly preferential orientation with molecular z along sample z axis
Exp.Ordering = +10; % strongly aligned such that molecular z axis is perpendicular to sample z axis

The built-in ordering function is very simple (see the pepper documentation for more information), but useful for many cases. Here is a simulation of a sample where the molecules are preferentially oriented such that the molecular z axis is parallel to the sample z axis: 

Sys.g = [2 2 2.2];
Sys.lwpp = 1;      % mT
Exp.mwFreq = 9.5;  % GHz
Exp.Ordering = +2;
pepper(Sys,Exp);

The second way to provide information about the orientational distribution is to write a separate function and supply to pepper as a function handle in Exp.Ordering. See the pepper documentation for details. The function takes as input the three Euler angles that transform the sample frame to the molecular frame (in analogy to Exp.MolFrame). Here is an example. Written as an anonymous function, the built-in orientational distribution is 

Exp.Ordering = @(alpha,beta,gamma) exp(lambda*plegendre(2,0,cos(beta)));

The EasySpin function plegendre returns the associated Legendre polynomial and lambda corresponds to the number given in Exp.Ordering.

Frequency-swept spectra

pepper, like the other cw EPR simulation functions garlic and chili, does field sweeps by default. However, you can ask them to simulate frequency-swept spectra as well.

For this, all you need to do is the following

Here is an example of a frequency-swept spectrum of a nitroxide radical: 

clear
Sys.g = [2.008 2.006 2.002];
Sys.Nucs = '14N';
Sys.A = [20 20 100];       % MHz
Sys.lwpp = 10;             % peak-to-peak line width, in MHz
Exp.Field = 340;           % static field, in mT
Exp.mwRange = [9.3 9.7];   % frequency range, in GHz
pepper(Sys,Exp);

By default, pepper returns the absorption spectrum (Exp.Harmonic=0) when you simulate a frequency-swept spectrum. To get the first or second derivative, change Exp.Harmonic to 1 or 2.

All other capabilities of pepper apply equally to frequency sweep and to field sweeps. For example, you can simulate crystals or multi-component spectra, and you can change the excitation mode. Importantly, just like field sweeps, frequency sweeps can be simulated using different methods: matrix diagonalization (Opt.Method='matrix'; highly accurate but potentially slow) or perturbation theory (Opt.Method='perturb', Opt.Method='perturb1', Opt.Method='perturb2'; less accurate but much faster).

Degeneracies

EasySpin is not able to handle systems with degenerate energy levels in an exact fashion when using matrix diagonalization. To make it possible to simulate EPR spectra of such systems, EasySpin automatically adds small random perturbations to non-zero Hamiltonian matrix elements. The level of this added noise can be adjusted in Opt.FuzzLevel. The default value is 1e-10 and should be adequate for most cases.

An alternative way to break degeneracies is to modify the spin system parameters. For instance, if there are two spins with identical g values, break the symmetry by adding a small perturbation to one of them: 

Sys.g = [2.01;  2.01];        % identical g values; Opt.FuzzLevel will break degeneracies
Sys.g = [2.01;  2.01+1e-7];   % slightly perturbed g values, to break degeneracies

If you break degeneracies by modifying spin system parameters, and in cases where you are looking for very small changes in the EPR line positions, it is advisable to turn off the addition of random shifts: 

Opt.FuzzLevel = 0;