chili

Simulation of field- and frequency-sweep cw EPR spectra of tumbling spin systems in the slow-motional regime.

Syntax
chili(Sys,Exp)
spec = chili(...)
[B,spec] = chili(...)
[nu,spec] = chili(...)
[B,spec,info] = chili(...)
[nu,spec,info] = chili(...)
... = chili(Sys,Exp,Opt)

See also the user guide on how to use chili.

Description

chili computes cw EPR spectra in the slow-motional regime. It uses rotational Brownian diffusion to model the rotational tumbling of the spin system. The simulation is based on solving the Stochastic Liouville equation in a basis of rotational eigenfunctions. chili supports arbitrary spin systems.

chili takes up to three input arguments

If no input argument is given, a short help summary is shown (same as when typing help chili).

Up to three output arguments are returned:

If no output argument is requested, chili plots the spectrum.

chili can simulate field-swept spectra as well as frequency-swept spectra. For field-swept spectra, specify Exp.mwFreq (in GHz) and Exp.Range (in mT); for frequency-swept spectra specify Exp.Field (in mT) and Exp.mwRange (in GHz).

Input: Spin system

Sys is a structure containing the parameters of the spin system. See the documentation on spin system structures for details. Most of the common spin system parameters are supported. However, chili does not support strains (gStrain, etc.).

For simulating a multi-component mixture, give a cell array of spin systems in Sys, 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.

In addition to the static spin system parameters, Sys should contain dynamic parameters relevant to the motional simulation. To specify the rate of rotational diffusion, give one of the fields tcorr, logtcorr, Diff, or logDiff.

tcorr
Rotational correlation time, in seconds.

For example,

Sys.tcorr = 1e-9;         % isotropic diffusion, 1 ns correlation time
Sys.tcorr = [5 1]*1e-9;   % axial anisotropic diffusion, 5 ns around x and y axes, 1 ns around z
Sys.tcorr = [5 4 1]*1e-9; % rhombic anisotropic diffusion

Instead of tcorr, Diff can be used, see below. The correlation time tcorr and the diffusion rate Diff are related by

tcorr = 1/(6*Diff);       % tcorr in s, Diff in rad^2 s^-1

logtcorr
Base-10 logarithm of tcorr, offering an alternative way to specify the rotational dynamics. Use this instead of tcorr for least-squares fitting with esfit.
Sys.logtcorr = -9;      % corresponds to tcorr = 1e-9 s
Diff
Principal values of the rotational diffusion tensor, in s-1 (equivalently, rad2 s-1).
logDiff
Base-10 logarithm of Diff. Use this instead of Diff for least-squares fitting with esfit.
DiffFrame
3-element vector [a b c] containing the Euler angles, in radians, describing the orientation of the rotational diffusion tensor frame relative to the molecular frame. DiffFrame gives the angles for the transformation of the molecular frame into the rotational diffusion tensor eigenframe. See frames for more details.

In addition to the rotational dynamics, convolutional line broadening can be included using Sys.lw or Sys.lwpp.

lwpp
1- or 2-element array of peak-to-peak (PP) line widths (all in mT for field-swept spectra, and in MHz for frequency-swept spectra).
lw
1- or 2-element array of full-width-half-maximum (FWHM) line widths (all in mT for field-swept spectra, and in MHz for frequency-swept spectra).

chili can model hindered rotation using an orientational potential, which can be specified in Sys.Potential.

Potential

An array of coefficients for the orientational potential, with four numbers per row: L, M, K, and λ. The potential is expressed as U(Ω) = - kB T ΣL,M,KλLM,KDLM,K, where the DLM,K are Wigner functions and λLM,K are the possibly complex-valued coefficients. L is a nonnegative integer, and both M and K are integers between -L and L.

Sys.Potential = [2 0 0 0.3];                  % L=2, M=K=0 coefficient is 0.3
Sys.Potential = [2 0 0 0.3; 2 0 2 -0.1];      % two coefficients

Since the overall potential must be real-valued, there is a symmetry relation between coefficients with the same L, |M| and |K|: λL-M,-K = (-1)M-KLM,K)*. Therefore, chili limits the input coefficients to: (1) those with positive K, (2) those with K=0 and positive M, (3) the one with M=K=0 which must be real-valued. All others are supplemented internally using the symmetry relation.

To plot the orientational potential, use the function oripotentialplot by calling oripotentialplot(Sys.Potential).

The frame of the orientational potential is assumed to be collinear with that of the rotational diffusion tensor (Sys.DiffFrame).

The orientational symmetry of the potential reflects the symmetry of the paramagnetic molecule and the symmetry of the nanoenvironment in which it is tumbling:

chili simulation performance varies greatly with the type of potential. The simulations are fastest if all potential terms have M=K=0, they are fast if all terms are M=0 or if all terms are K=0, and they are very slow if any term has both M and K different from zero.

For concentrated solutions, it is possible to include Heisenberg exchange:

Exchange

Effective Heisenberg spin exchange rate, in inverse microseconds. Implements a simple contact-exchange model, see eq. (A27) from Meirovitch et al, J.Chem.Phys.77, 3915-3938. See also Freed, in: Spin Labeling (ed. L.J. Berliner), 1976, p.68.

Input: Experimental parameters

The experiment structure Exp contains all parameters relating to the experiment. These settings are identical for all cw EPR simulation functions (pepper, chili, garlic). See the page on cw EPR experimental parameters.

Additionally, the following field is supported:

Ordering scalar (default: zero) or function handle

If a number is given in this field, it specifies the orientational distribution of the supporting biomacromolecule in the sample. If not given or set to zero, the distribution is isotropic, i.e. all orientations occur with the same probability.

If a number is given, the orientational distribution is non-isotropic and computed according to the formula P(β) = exp(-U(β)) with U(β) = -λ(3 cos2β - 1)/2, where β is the angle between the sample z axis (z axis of the sample frame) and the molecular z axis ( axis of the molecular frame), and λ is the number specified in Exp.Ordering.

Typical values for λ are between about -10 and +10. For negative values, the orientational distribution function P(β) is maximal at β = 90° (preferential orientation of molecular z axis in the sample xy plane), for positive values it is maximal at β = 0° and β = 180° (preferential alignment of the molecular z axis with the sample z axis). The larger the magnitude of λ, the sharper the distributions.

To plot the orientational distribution for a given value of λ, use

lambda = 2;
beta = linspace(0,pi,1001);
U = -lambda*plegendre(2,0,cos(beta));
%U = -lambda*(3*cos(beta).^2-1)/2;  % equivalent
P = exp(-U);
plot(beta*180/pi,P);

If Exp.Ordering is a function handle, pepper calls the function to obtain the orientational distribution. The function must accept either one (beta) or two input arguments (beta and gamma), in radians. The function must accept vector arguments and return a vector P containing probabilities for each orientation, that is P(k) is the probability of finding the spin centers with orientation (relative to the sample frame) specified by beta(k) and gamma(k). Here is an example with an anonymous function:

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

If the function cannot be written as a singe line, define a separate function and provide it via its function handle. For example, here is a function oridist defining a von Mises-Fisher distribution

function P = oridist(beta,gamma)
gamma0 = 0;
beta0 = pi/4;
kappa = 10;
mu = ang2vec(gamma0,beta0);
x = ang2vec(gamma(:),beta(:));
P = kappa/(4*pi)/sinh(kappa)*exp(kappa*mu'*x);
P = reshape(P,size(beta));
end

Include this in the simulation via Exp.Ordering = @oridist.

When providing a function for Exp.Ordering, make sure that the symmetry used in the simulation is the same or lower than the symmetry of the distribution. Otherwise, incorrect spectra are obtained. In cases of doubt, set Opt.GridSymmetry='C1'. This always gives correct results.

If the orientational distribution provided via Exp.Ordering is very narrow, increase the number of knots in Opt.GridSize.

Input: Simulation options

Opt, the options structure, collects all settings relating to the algorithm used and the behavior of the function. The most important settings are LLMK, highField, pImax, and GridSize, since they determine the basis size and the number of orientations for powder averaging. If any of these settings is too small, the spectrum will not be converged.

LLMK
4-element vector [evenLmax oddLmax Mmax Kmax]
Specifies the rotational basis size by giving the maximum values for, in that order, even L, odd L, M, and K. All four numbers must be non-negative. The maximum values for M and K must be less than or equal to the maximum value of L. The larger the numbers, the larger the rotational basis.
If this field is not specified, chili uses a default basis. This is adequate only for common X-band spectra of nitroxides. In general, the basis needs to be larger for slower motions and can be smaller for faster motions. It is strongly advised to increase these four numbers until the simulated spectrum is converged. One approach is to start from [2 0 0 0] and increase each of the four numbers in turn, run the simulation, and see if the spectrum changes. Once the spectrum doesn't change anymore, the basis is sufficiently large.
highField
-1 (default) or +1
Determines which electron spin transitions to include in the spin basis. If set to false, includes all transitions with pS = mS'-mS'' = -1 and +1. This is applicable in all cases. In certain situations, one can get away with setting highField to true. This indicates the high-field approximation, and that only transitions with pS = +1 are included. This reduces the basis and speeds up simulations, but can result in incorrect spectra under certain circumstances. It is important to always check such simulations against highField=false.
pImax
vector, one element per nucleus
Determines for each nucleus the maximum nuclear coherence order |ΔmI| to include in the spin basis. 0 corresponds to only allowed EPR transitions, 1 includes additionally forbidden EPR transitions with |ΔmI|=1, etc. By default, all possible orders are included. By reducing pImax, the spin basis is truncated, and the calculation runs faster. For nitroxides, truncation to Opt.pImax=1 often works and can yield spectra fairly close to the correct one. Truncation to Opt.pImax=0 rarely works and is generally not advisable.
pImaxall
Gives the maximum total nuclear coherence order (summed over all nuclei) to include in the spin basis. This truncation is applied in addition to the one given in Opt.pImax.
GridSize
Size of orientational grid used in a powder simulation (i.e. if the protein or other director does not have a single orientation, but is orientationally disordered). Unless Exp.SampleFrame is set, a powder simulation is performed. The default value for GridSize is 19. Increase this value if the orientational potential coefficients in Sys.Potential are large.
GridSymmetry
'Dinfh' (default), 'D2h', 'C2h' or 'Ci'
Determines the symmetry of the orientational grid for powder simulations. '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). The default 'Dinfh' is correct only for high-symmetry spin systems in high-symmetry potentials. In general, it is advised to run a simulation with Ci, which is slower, but correct for any situation.
LiouvMethod
This specifies which method is use to construct the Liouville matrix (i.e. the Hamiltonian and the relaxation superoperator). The two possible values are 'fast' and 'general'. The fast method is very fast, but limited to one electron spin with S=1/2 coupled to up to two nuclei and to orientational potentials with even L≤4, zero M, even K ≤2, and real coefficients. On the other hand, the general method works for any spin system and any form of potential, but is significantly slower. By default, chili uses the fast method if applicable and falls back to the general method otherwise.
FieldSweepMethod
Specifies how to calculate the field sweep. There are three possible settings: 'explicit', 'approxinv', and 'approxlin'. 'explicit` indicates that the spectral intensity should be calculated separately and explicitly for each field value. This is the accurate method. In many cases, using 'approxlin' or 'approxinv' results in much faster simulations. For these settings, a frequency-swept spectrum is calculated first and then converted to the field domain. This procedure is only approximate. 'approxinv' works best if the g anisotropy dominates the spectrum, in all other cases (e.g. dominant hyperfine or zero-field splittings), 'approxlin' is preferable.
separate '' (default), 'components'

Determines whether to return the total spectrum or a list of subspectra. If set to '', the total spectrum is returned. If set to 'components', spec is a matrix with the subspectra of all components (including isotopologues). Each row in spec is one subspectrum.

PostConvNucs
This specifies which nuclei should be excluded from the Stochastic Liouville equation (SLE) simulation and only included in the final spectrum perturbationally, via post-convolution of the SLE-simulated spectrum with an isotropic stick spectrum of the nuclei marked for post-convolution. E.g. If Sys.Nucs = '14N,1H,1H,1H' and Opt.PostConvNucs = [2 3 4], the only the nitrogen is used in the SLE simulation, and all the protons are included via post-convolution.

Post-convolution is useful for including the effect of nuclei with small hyperfine couplings in spin systems with many nuclei that are too large to be handled by the SLE solver. Nuclei with large hyperfine couplings should never be treated via post-convolution. Only nuclei should be treated by post-convolution for which the hyperfine couplings (and anisotropies) are small enough to put them in the fast-motion regime, close to the isotropic limit, for the given rotational correlation time in Sys.tcorr etc.

Solver
'L', '\', 'E', 'B', or 'C'
Optionally, specifies the linear solver to use. Possible values are 'L' (Lanczos tridiagonalization with continued-fraction expansion), '\' (MATLAB's backslash solver), 'E' (via eigenvalues), 'B' (biconjugate gradient stabilized method), and 'C' (conjugate gradient tridiagonalization with continued-fraction expansion). If not given, a solver is selected automatically.
Verbosity
0 (default), 1
Determines how much information chili prints to the command window. If Opt.Verbosity=0, is is completely silent. 1 prints details about the progress of the computation.
Example

The cw EPR spectrum of a nitroxide radical tumbling in solution on a time scale of nanoseconds can be simulated with the following lines.

clear
Sys.g = [2.008 2.0061 2.0027];
Sys.Nucs = '14N';
Sys.A = [16 16 86];   % MHz
Sys.tcorr = 3e-9;     % 3 ns
Exp.mwFreq = 9.5;     % GHz
chili(Sys,Exp);
Algorithm

chili solves the Stochastic Liouville equation (SLE). It represents the orientational distribution of the spin system using normalized Wigner rotation functions DLK,M(Ω) with -L≤K,M≤L as basis functions. The number of basis functions is determined by maximum values of even L, odd L, K and M. The larger these values, the larger the basis and the more accurate the spectrum.

chili computes EPR line positions to first order, which is appropriate for most organic radicals. It is inaccurate for transition metal complexes, e.g. Cu2+ or VO2+. For the diffusion, both secular and nonsecular terms are included.

If the spin system has S=1/2 and contains no more than two nuclei, chili by default uses a fast method to construct the Liouvillian matrix that uses code based on the software from the Freed lab at Cornell (Opt.Method='fast'). For all other cases, a general code is used to construct the Liouvillian matrix.

Post-convolution works as follows: First, the SLE is used to simulate the spectrum of all nuclei except those marked for post-convolution. Next, the isotropic stick spectrum due to all post-convolution nuclei is simulated and convolved with the SLE-simulated spectrum to give the final spectrum.

For full details of the various algorithms see

See also

esfit, fastmotion, garlic