EPR Simulations with EasySpin
Installing Software
EasySpin is a MATLAB toolkit, so you must first install MATLAB. You can sign up for a trial version by entering your email here .
1) Create your MathWorks account and follow the prompts...
2) Once you're on the download page, choose the installer that meets your OS.
3) Once this is downloaded, click on it to run the installation software.
4) Once prompted, sign into your MathWorks account and select the license you would like to download and install for the trial.
5) Select the products you would like to download and install. (I left all of them selected.)
6) Let the installation complete... You now have MATLAB.
To get EasySpin, you must download it here (click on the "Get Latest Version" button).
7) Now, "Extract" or Unzip the files you have downloaded. (Location doesn't matter, entirely up to you.)
8) Now, launch MATLAB. You have to "connect" EasySpin to MATLAB...
9) To do this, you must click on Home->Set Path (see picture below)
10) Click on "Add Folder" and select the EasySpin subfolder, called easyspin (easyspin-x.y.z/easyspin). Then click Save.
This should have made it possible for MATLAB to communicate with EasySpin. You can test this by typing easyspin into the MATLAB command prompt... Happy simulating!
Overview
This page serves as a resource to assist those trying to simulate EPR spectra with the EasySpin MATLAB package.
We should start out with a simple example, to see what kind of input EasySpin requires. Then, we will talk about the details regarding the parameters encountered in the example.
A simple example simulation:
Sys = struct('g',2,'Nucs','1H','n',2,'A',15,'lwpp',[0 0.05]);
Exp = struct('mwFreq',9.5);
garlic(Sys,Exp);
This should yield something that looks like:
Importing Data
EasySpin has a function for doing this, called eprload(). The method used to import your data depends on the format in which your experimental data is saved:
- For Modern Bruker Spectrometer Data:
[B,spc] = eprload('mydata');
- For Older Bruker Spectrometer Data:
[B,spc] = eprload('myolddata.spc');
[B,spc] = eprload('myolddata.par');
The function changes when you are trying to load a text file of data:
- For Data Stored in Text Files:
[B,spc] = textread('mydata.txt','%f %f');
- Text Files with Header Information:
[B,spc] = textread('mydata.txt','%f %f','headerlines',17)
Creating a Simulation
The central concept to constructing a simulation with EasySpin is defining a Spin System (Sys) and Experimental Parameters (Exp).
Spin System
There are a variety of components that factor into a Spin System.
Electron and Nuclear Spins
The spins can be defined with Sys.S() according to the following:
Sys.S = 1/2; % one electron spin with S=1/2
Sys.S = 5/2; % an S=5/2 spin
Sys.S = [1/2, 1/2]; % two S=1/2 spins
Sys.S = [1, 1, 1/2]; % two S=1 spins and one S=1/2 spin
If this is not specified, the default value is S = 1/2.
The Nuclei can be specified using Sys.Nucs() according to the following:
Sys.Nucs = '1H'; % one hydrogen
Sys.Nucs = '63Cu'; % a 63Cu nucleus
Sys.Nucs = '59Co,14N,14N'; % a 59Co and two 14N nuclei
You can specify an arbitrary number of nuclei for this system. In order to specify the number of each type of nucleus, you must use Sys.n() :
Sys.Nucs = '1H,13C'; % one class of 1H and one class of 13C
Sys.n = [2 3]; % 2 protons and 3 carbon-13 spins
G-value
As you're on this page, you are probably aware that there are some considerations to be made about the kind of system (isotropic vs. anisotropic) we are simulating when specifying the g-value. For all purposes here, we will assume that it will always be isotropic. To specify this you must use Sys.g():
Sys.g = 2.0023;
This actually is not necessary, since the default is set to 2.0023... for all inputs. Still worth mentioning!
Hyperfine Coupling Constants and Number of Equivalent Nuclei
These are specified in MHz in EasySpin. No worries, though because EasySpin has a function that allows us to convert mT (or Gauss) into MHz:
A_MHz = mt2mhz(A_mT); % mT -> MHz conversion
A_MHz = mt2mhz(A_G/10); % G -> MHz conversion (1 G = 0.1 mT)
Consider the following system, in which you have a hydrogen atom with a 10 MHz coupling to the unpaired electron and a 13C atom with a 12 MHz coupling:
Sys.Nucs = '1H,13C'; % Defines the Hydrogen nucleus and the Carbon-13 nucleus
Sys.A = [10 12]; % MHz... Defines the HFCC
To specify the number of equivalent nuclei, Sys.n() is used, as shown in the following examples:
For a single proton:
Sys.Nucs = '1H';
Sys.n = 1;
Sys.A = 5.3;
For two equivalent protons (with identical hyperfine coupling constants):
Sys.Nucs = '1H';
Sys.n = 2;
Sys.A = 5.3;
For something more complex (such as the naphthalene radical anion, which contains two sets of 4 equivalent nuclei):
Sys.Nucs = '1H,1H';
Sys.n = [4 4];
Sys.A = [-14 -5];
Line Width
To specify the line width, EasySpin uses Sys.lwpp(). This gives the peak-to-peak line width. It should be noted here that the unit used for Line Widths is mT!!! For peak-to-peak line widths:
Sys.lwpp = [0 0.03]; % Gaussian and Lorentzian peak-to-peak line width, mT
There is another option Sys.lw():
Sys.lw = [0 0.05]; % Gaussian (1st input) and Lorentzian (2nd input) FWHM (full width at half height), mT
Experiment Parameters
There are a variety of fields that you can/must specify before you can get your simulated spectrum.
Microwave Frequency
The microwave frequency must be given in units of GHz. This can easily be defined with the following:
Exp.mwFreq = 9.385; % X-band
Field Range/Sweep Width
This is an optional step. If omitted, EasySpin will attempt to determine a Sweep Width large enough to accommodate the entire spectrum, but this can fail.
There are two ways to enter the magnetic field sweep width. Either give the center field and the sweep width (in mT) in Exp.CenterSweep():
Exp.CenterSweep = [340 80]; % in mT
As an alternative, you can specify the lower and upper limit of the sweep range (again in mT) in Exp.Range():
Exp.Range = [300 380]; % in mT
Number of Points
Again, an optional step. It is not required to be a power of 2 for EasySpin.
Exp.nPoints = 5001;
Harmonic
You can set the harmonic that EasySpin will calculate, meaning either the Absorption or a derivative. This is specified with Exp.Harmonic():
Exp.Harmonic = 0; % absorption spectrum, direct detection
Exp.Harmonic = 1; % first harmonic (default)
Running the Simulation
Once you are ready to run the simulation, you simply have to call the garlic() function on your Spin System variable (Sys) and Experimental Set-up Variable (Exp):
garlic(Sys,Exp);
Multiple Species in one Simulation
EasySpin can handle this fairly well, by simply defining a second Spin System:
Sys1.g = 2;
Sys1.lwpp = 1;
Sys1.weight = 2;
Sys2.g = 2.1;
Sys2.lwpp = 0.8;
Sys2.weight = 1;
Exp.mwFreq = 9.5;
Exp.Range = [300 360];
garlic({Sys1,Sys2},Exp);
Fitting EPR Spectra with EasySpin
The Fitting Function
In order to fit an EPR Spectrum, EasySpin uses a function called esfit(). The general structure of this command is as follows:
esfit('simulation function',spc,Sys0,Vary,Exp)
An explanation of the parameters of the esfit() function:
- 'simulation function' is the EasySpin simulation function that will be used for the fit. For fast-motion, isotropic EPR spectra the function will be garlic as discussed above.
- spc is the array containing the experimental spectral data
- Sys0 is a structure containing the magnetic parameters of the spin system. These parameters are used as first guesses in the fitting procedure. For multiple systems, this will be a list: {Sys1, Sys2}.
- Vary is a structure similar to Sys0 containing the search ranges of the parameters that should be fitted. For a multicomponent fit, this is again a list of structures, one for each component: {Vary1, Vary2}. This should have the same defined parameters as Sys0.
- Exp is a structure containing the experimental parameters for the simulation function, as discussed above.
Loading the Experimental Spectrum
The first step is to import the experimental data. This is done in the same fashion as we did earlier.
[B, spc] = textread('filename.txt', '%f %f', 'header', number of header lines)
Initial Set of Parameters
Sys0.g = 2.0; Sys0.lw = 1;
Parameter Search Range
This allows us to tell esfit what to change in Sys0 by utilizing the Vary class.
Vary.lw = 0.3
These variations are the original value defined in Sys0 +/- the value defined in Vary.
You should vary no more than 4 parameters at once. As the number of degrees of freedom increases, the difficulty of the fit dramatically increases.
Launching the Fitting Graphical User Interface (GUI)
Now that you have defined all of the parameters needed to run the fit of your experimental spectrum, you simply need to call the esfit() function:
esfit('simulation function',spc,Sys0,Vary,Exp)
This will launch the GUI, pictured below. This can be very daunting at first, but we will take a look at what all of these parameters mean.
Fitting Methods
Here, we will take a look at what the various optimization algorithms and target functions do.
Optimization Algorithms
There are five optimization algorithms in EasySpin:
- Neider/Mead Downhill Simplex
- Levenberg/Marquardt Algorithm
- Monte Carlo Random Search
- Genetic Algorithm
- Particle Swarm/Systematic Grid Search
Neider/Mead Downhill Simplex and Levenberg/Marquart Algorithm are both local search algorithms, which begin with the given starting parameters and then work their way down a nearby valley of the parameter space to find a minimum. Both are relatively fast methods, but there are some general differences. The downhill simplex is slower than the Levenberg/Marquardt, but it is less likely to get stuck in a local minimum.
The remaining three algorithms are global search methods: they do not have a starting set of parameters, but use a variety of sets distributed over the entire parameter search space.
The Monte Carlo method simply performs a series of random trial simulations and then picks the best one. It is very inefficient. The systematic grid searches better because it covers the parameter space with a grid and then does simulations for each knot of the grid in random order. No point is simulated twice and this is a more efficient method than the Monte Carlo search. But, if the minimum is between two grid points, it will never be found.
The genetic algorithm makes simulations for several (N) parameter sets (called a population), computes the fitting error (called the fitness) for all of them, and then proceeds to generate N new parameter sets from the old ones using mechanisms like mutation, cross-over, and reproduction. This way, a new generation of parameter sets is procreated, as in biological evolution. The benefit of this algorithm is that if a good parameter is encountered, then it is likely to propagate down the generations and across the population.
Target Function
The second important setting is the target function. 'esfit' computes the error of the simulated spectrum using the root-mean-square-deviation (rmsd, which is the square root of the mean of the squared standard deviations), where the deviation is the difference between the simulated and experimental spectra.
Fitting speed can often be increased if the spectra are not directly used, that is, if their integrals or similar transforms are used instead. The options in EasySpin are as follows
- Data As Is
- Integral
- Double Integral
- First Derivative
- Fast Fourier Transform
A Guided Tutorial
In order to see how all of the previous topics are involved in the EPR Simulation Process, we will walk through a guided tutorial where we will simulate real experimental data.
This tutorial requires the following input files:
General Structure of the Template File
close all; This line closes any pre-existing MATLAB Figures.
clear all; This line clears all pre-existing variable definitions.
[B,spec1] = textread('filename.txt', '%f %f');
This line defines the magnetic field (B) in units of the experimental spectrum (spec1). Recall that EasySpin uses units of mT, so you will have to convert your experimental spectrum if it was saved in units of Gauss. Here %f is referred to as a string formatter, allowing the textread() function to know that the experimental data is in "floating point decimal" format - indicating that the data is numerical with a variable position of the decimal point.
The next lines define our simulation system - here identified as M :
M.S = System Spin; Define the spin of our system.
M.Nucs = 'nuclei'; Define any nuclei present/coupled to the spin.
M.g=[g-value(s)]; Define the principle g value or values for the system. (It's not always 2.0023.....)
M.A =[Hyperfine Coupling Constant (in mHz)]; Define any hyperfine couplings present in the system (should have the same shape as M.Nucs)
M.lwpp=[0.1 0.1]; % Defaults to units of what we are sweeping (mT/MHz field/frequency) Define our line width broadening -
% Note that this is in the format of a Voigt curve with [Gaussian_Component Lorentzian_Component].
Next we have to set-up the Experimental Parameters for the simulation:
Exp=struct('mwFreq',x.xx,'nPoints',xxxx,'Range',[xx.xx xx.xx]); Define our magnetic field window for the simulation - which can be easily obtained from the experimental set-up.
Note: nPoints should match the number of points in the experimental data file.
Exp.Temperature = x; Define the temperature at which the experimental spectrum was collected.
We then run the simulation using a specific regime function (garlic, pepper, etc.):
[B,simul1]=regimefunction(M,Exp);
Finally, we can apply a baseline correction to our experimental spectrum or an intensity correction to our simulated spectrum:
y=spec1+b % Baseline Correction/y-shift
y2=simul1*k % Intensity correction for simulated data.
Important Note: the template MATLAB file provided above is currently saved as a '.txt' file (so I could upload it to the Wiki server) - you will need to change this extension to '.m' in order to run in MATLAB.
Simple Example: DPPH
DPPH is a common standard used in EPR. This has a well-characterized spectrum, which makes identification of the principle g-value and hyperfine coupling constant very easy.
We will now walk through the steps of simulating this spectrum.
- Create a folder named "DPPH". Add to this folder both the experimental data file (DPPH.txt) and a copy of the MATLAB template file.
- Open your template file and save it as something like DPPH_simulation.m
- We will now adjust some lines of the original template to begin our simulation.
- Change filename.txt to DPPH.txt
[B,spec1] = textread('DPPH.txt', '%f %f');
- We are dealing with a lone free radical in DPPH without hyperfine coupling, so we must change System Spin to 1/2 and comment out M.Nucs and M.A
M.S = 1/2;
%M.Nucs = 'nuclei';
%M.A =[Hyperfine Coupling Constant (in mHz)];
- DPPH has a well-characterized principle g-value of 2.0036, so insert this value
M.g=[2.0036];
- Play with the peak-peak line broadening.
M.lwpp=[0.1 0.1];
- Set the experimental parameters to:
Exp=struct('mw Freq',x.xx,'nPoints',xxxx,'Range',[xx.xx xx.xx]);
- Set the experimental temperature to 298 K:
Exp.Temperature = 298;
- We are in the fast-motion, liquid regime, so we will use garlic
[B,simul1]=garlic(M,Exp);
- Play with the baseline and intensity corrections. (Hint: This spectrum requires no baseline and has an intensity correction of 0.05.)
Your final output should look like the following figure:
Example 2: TEMPO
TEMPO is a nitroxide that is commonly encountered in the field of EPR. Its spectrum serves as an excellent example of an isotropic signal displaying hyperfine coupling. This adds another layer of difficulty to the procedure used in the DPPH spectrum above. The hyperfine coupling arises from delocalization of the radical across both the oxygen and nitrogen atoms. Nitrogen-14 has nuclear spin of I = 1, which accounts for the three lines in the spectrum. Note that there are also some satellite peaks located on each side of each of the three peaks.
Now, let's get started with the simulation:
- Create a folder named "TEMPO". Add to this folder both the experimental data file (TEMPO.txt) and a copy of the MATLAB template file.
- Open your template file and save it as something like DPPH_simulation.m
- Modify the code and see how well you can simulate the experimental spectrum, taking note of the following experimental parameters:
- Microwave Frequency: 9.869 GHz
- Number of Points: 485
- Field Range: 347 - 356.5 mT
- Temperature: 298 K
- TEMPO is an organic radical, so you can safely guess that the principle g-value will be around 2.0. If you'd rather handle this directly, you can utilize the following relationship between the magnetic field, microwave radiation, and g-value to calculate a g-value to use for your guess:
- As previously mentioned, TEMPO exhibits hyperfine coupling. This means that you will have to define a nucleus in the M.A line. Recall from above that EasySpin thinks about Hyperfine Couplings in units of MHz. If you prefer to think about these values in terms of Magnetic Field units (Gauss or milliTesla), you can use the mt2mhz() to wrap your input as follows:
M.A = [mt2mhz(A in mT)];
Note that you can estimate the hyperfine splitting directly from the spectrum and then convert to MHz with the following equation:
- Play with the Hyperfine Coupling Constant (A), line width (lwpp), g, and the scaling factor until you have a sufficient fit of the main triplet signal. Once you've successfully simulate the triplet, you may move onto the satellite peaks.