This user guide explains how to simulate trajectory-based cw EPR spectra of using EasySpin's function cardamom. It is assumed that you are familiar with the basics of MATLAB, esp. with structures.

• Running the simulation
• Specifying the static parameters
• Dynamic parameters and model
• The orientational potential
• Basic experimental settings
• Time axis and grid parameters
• Powder vs single orientation vs global diffusion
• Simulation options
• MD settings

Trajectory-based cw EPR spectra are computed by the EasySpin function cardamom.

cardamom(Sys,Exp,Par)

It is called with three arguments. The first argument Sys tells cardamom all about the static and dynamic parameters of the spin system. The second argument Exp gives the experimental parameters. The third argument Par gives the simulation parameters.

If no output argument is given, cardamom plots the computed spectrum. But it can also return one, two, or four outputs. (Don't forget the semicolon at the end of the line to suppress output to the command window.)

[Field,Spec] = cardamom(Sys,Exp,Par);
[Field,Spec,FID,time] = cardamom(Sys,Exp,Par);

The outputs Field and Spec are vectors containing the magnetic field axis and the spectrum, respectively, whereas FID and time are vectors containing the FID signal and the time axis. If these are requested, cardamom does not plot the spectrum.

Doing a simulation only requires a few lines of code. A simple example is

Sys = struct('g',[2.008,2.006,2.003],'Nucs','14N','A',[20,20,85]);
Sys.tcorr = 4e-8;
Exp = struct('mwFreq',9.5);
Par = struct('dtSpin',2e-9,'dtSpatial',2e-9,'nSteps',100,'nTraj'400);
cardamom(Sys,Exp,Par);

The first line defines the spin system, a nitroxide radical with anisotropic g and A tensors. The second line gives the rotational correlation time of the radical. The third line specifies experimental parameters, here only the microwave frequency is needed (the magnetic field range is chosen automatically). The fourth line defines the simulation parameters and trajectory properties. The fifth line calls the simulation function, which plots the result. Copy and paste the code above to your MATLAB command window to see the graph.

Of course, the names of the input and output variables don't have to be Sys, Par, Exp, Field, and Spec. You can give them any name you like, as long as it is a valid MATLAB variable name, e.g., FremySaltSolution or QbandExperiment. For convenience, thoughout this tutorial, we will use short names like Sys and Exp.

The first input argument specifies the static parameters of the paramagnetic molecule. It is a MATLAB structure with various fields giving values for the spin system parameters.

Sys.g = [2.008,2.006,2.003];
Sys.Nucs = '14N';
Sys.A = [20,20,80];  % MHz

The first line defines the g tensor of the spin system via its three principal values. cardamom always assumes a single unpaired electron spin S=1/2 if no spin is given.

The field Sys.Nucs contains a character array giving all the magnetic nuclei in the spin system, a nitrogen-14 in the above example. If Sys.Nucs is not specified, cardamom always assumes a ¹⁴N nitroxide spin label. If a simulation with a different type of nucleus is desired, then a different simulation method called 'ISTOs'needs to be declared in Opt.Method. For more details, please see cardamom.

Sys.A gives the hyperfine couplings in MHz (megahertz), with three principal values on a row for each of the nuclei listed in Sys.Nucs.

Remember that cardamom (and other EasySpin functions, too), take the hyperfine coupling values to be in MHz. Often, values for hyperfine couplings are given in G (Gauss) or mT (millitesla), so you have to convert these values. For g = 2.00232, 1 G corresponds to 2.8025 MHz, and 1 mT corresponds to 28.025 MHz. The simplest way to convert coupling constants from magnetic field units to MHz is to use the EasySpin function unitconvert:

A_MHz = unitconvert(A_mT,'mT->MHz');    % mT -> MHz conversion
A_MHz = unitconvert(A_G/10,'mT->MHz');  %  G -> MHz conversion (1 G = 0.1 mT)

The spin system structure also collects parameters relating to the dynamics of the paramagnetic molecules.

There are two types of local dynamics models in cardamom that can be used to simulate a cw EPR spectrum: stochastic, which simulates rotational diffusion, and jump, which simulates jumps between a set of discrete orientations. Each of these models can be declared using the field Par.Model.

For a rotational diffusion model (Par.Model='diffusion'), there are several alternative ways to specify the rate of isotropic rotational diffusion: either by specifying tcorr, the rotational correlation time in seconds

Sys.tcorr = 1e-7;   % 10^-7 s = 100 ns

or by giving its base-10 logarithm

Sys.logtcorr = -7;   % 10^-7 s = 100 ns

or by specifying the principal value of the rotational diffusion tensor (in s^-1)

Sys.Diff = 1e9;  % 1e9 s^-1 = 1 ns^-1

or by giving its base-10 logarithm

Sys.logDiff = 9;   % 1e9 s^-1 = 1 ns^-1

Diff and tcorr are related by

Diff = 1/6/tcorr;

The input field Diff can be used to specify an axial rotational diffusion tensor, by giving a 2-element vector with first the perpendicular and the the parallel principal value:

Sys.Diff = [1 2]*1e8;  % in hertz

For a jump model (Par.Model='jump'), the field TransRates can be used to specify the transition rate matrix (for two states in this case):

Sys.TransRates = [-0.5,  0.2;
                   0.5, -0.2]*1e8;  % in hertz

Alternatively, the transition probability matrix TransProb can be given for a stochastic jump model:

Sys.TransRates = [0.3, 0.36;
                  0.7, 0.64];

Note that Orientations, a set of Euler angles whose size is equal to the number of states (number of rows or columns in either of the above matrices), must also be given for a jump model:

Sys.Orientations = [   0,    0, 0;
                    pi/4, pi/2, 0];

The fields Sys.lwpp and Sys.lw has the same meaning as for the other EPR simulation functions (garlic, chili, pepper). You can use them to specify an additional convolution Gaussian and/or Lorentzian broadening in mT (Sys.lwpp for peak-to-peak, and Sys.lw for full width at half maximum).

Sys.lwpp = [0.05 0.01];   % [GaussianPP, LorentzianPP] in mT

For the reliability of the simulation algorithm we recommend to add a small residual Lorentzian broadening in Sys.lwpp.

cardamom can also include an orientational potential in the simulation. It is specified in the field Potential in the spin system structure. A Wigner D-function expansion can be used to specify the orientational potential by giving the L, M, and K values in the first three columns and the coefficient in the fourth column. Each term in the expansion then corresponds to each row in the matrix. For example:

Nx.Potential = [2, 0, 0,  0.3;
                2, 0, 2, -0.2];

indicates λ²_0,0 = 0.3 in the first row and λ²_0,2 = -0.2 in the second row.

Note that an orientational potential can only be used for a stochastic rotational diffusion model, not a jump model.

The second input argument, Exp, collects all experimental 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 enter 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

On many cw EPR spectrometers, the field range is specified using center field and sweep width, so Exp.CenterSweep is often the more natural choice.

Exp.CenterSweep and Exp.Range are only optional. If both are omitted, EasySpin tries to determine a field range large enough to accommodate the full spectrum. This automatic ranging works for most common systems, but fails in some complicated situations. EasySpin will issue an error when it fails.

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

Exp.nPoints = 5001;

You can set any value, unlike some EPR spectrometers, where often only powers of 2 are available (1024, 2048, 4096, 8192).

Since cardamom calculates cw EPR spectra by first simulating spin dynamics in the time domain, one must specify the parameters associated with time integration. This includes the spatial dynamics time step Par.dtSpatial associated with molecular motion, e.g. rotational diffusion of a spin label:

Par.dtStoch = 1e-9; % spatial propagation

Additionally, the spin dynamics time step Par.dtSpin and number of steps Par.nSteps are required to simulate the FID:

Par.dtSpin = 1e-9;  % spin propagation
Par.nSteps = 300;

Finally, it is necessary to specify the dynamics model type for the simulation (see the technical documentation for the full range of options):

Par.Model = 'diffusion'; % simulate stochastic rotational diffusion

Regarding orientational grid information, cardamom will try to determine the best parameters to simulate a cw EPR spectrum. However, depending on the desired model, sometimes it is necessary tell cardamom how the orientational grids should be handled.

In a frozen solution of spin-labeled protein, the protein molecules assume all possible orientations. For slow-motion spectra, the situation is more complex. The spin label will undergo "local" rotational diffusion, which may or may not be under the influence of an orientational potential. The protein itself, depending on its size and environment, might also be undergoing "global" rotational diffusion, which would couple to the local diffusion of the spin label. Otherwise, the protein might be "fixed" on the timescale of a cw EPR measurement. For the latter case, if the spin label is under the influence of an orientational potential, the distribution of protein orientations must be taken into account. If there is no potential, then it is sufficient to compute only one orientation, as the spectra from all orientations are identical.

The summation of slow-motion spectra over all possible orientations of an immobile protein ("director") is historically called the MOMD (microscopic order macroscopic disorder) model. This is equivalent to a powder spectrum. In cardamom, the powder spectrum is simulated whenever you specify an orientational potential. To get a single-crystal spectrum, i.e. the slow-motion spectrum for the nitroxide attached to a rigid protein with a single orientation, give the crystal orientation in Par.Orients. Par.Orients contains the Euler tilt angles (in radians), between the lab frame (which is lab-fixed) and the frame of the orienting potential (which is fixed to the protein).

When cardamom performs a powder simulation, it takes the number of orientations to include from Par.nOrients. Often, Par.nOrients does not have to be changed from its default setting, but if the spectrum does not appear to be smooth, Par.nOrients can be increased. Of course, this also increases the simulation time. For quick simulations, Par.nOrients should be minimized.

The third input structure, Opt, collects parameters related to the simulation algorithm.

The field Verbosity specifies whether cardamom should print information about its computation into the MATLAB command window. By default, its value is set to 0, so that cardamom is silent. It can be switched on by giving

Opt.Verbosity = 1;     % print information

Since cardamom simulates cw EPR spectra based on time-domain information, molecular dynamics (MD) trajectories of spin-labeled proteins can also be used to simulate spectra. To use this feature, the MD simulation output first must be processed to extract the information that cardamom needs to perform the simulation.

EasySpin can interface with MD simulation output using the function mdload. To do this, mdload needs to know about the trajectory output files using TrajFile and the relevant atomic information using AtomInfo. The output is then stored in the structure variable MD:

MD = mdload(TrajFile, AtomInfo);

TrajFile can be a single trajectory output file or a list of files. AtomInfo tells mdload how to look for spin label and protein information in TrajFile. For more details on how to do this, please see the technical documentation mdload.

Par.Model = 'MD-direct';
[Field,Spec] = cardamom(Sys,Exp,Par,Opt,MD);

By default, cardamom uses the orientations directly from the MD simulation trajectories to simulate the spectrum. For other trajectory usage options, see the technical documentation page for cardamom.

Note that EasySpin currently supports only NAMD and CHARMM formats (.PSF for topology files and .dcd for trajectory files). Please let the developers know if there is another file format that you would like to be supported.