Simulation of field- and frequency-sweep cw EPR spectra of tumbling spin systems in the slow-motional regime.
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
.
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
Sys
: static and dynamic parameters of the spin system
Exp
: experimental parameters
Opt
: options and settings
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:
B
or nu
: magnetic field axis, in mT, or frequency axis, in GHz
spc
: calculated spectrum
info
: structure containing details about the calculation
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).
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
txyz
: isotopic diffusion, corresponding to [txyz txyz txyz]
[txy tz]
: anisotropic diffusion with axial diffusion tensor, corresponding to [txy txy tz]
[tx ty tz]
: anisotropic diffusion with rhombic diffusion tensor
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
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
Dxyz
: isotopic diffusion tensor, corresponding to [Dxyz Dxyz Dxyz]
[Dxy Dzz]
: axial tensor, corresponding to [Dxy Dxy Dzz]
[Dxy Dxy Dzz]
: rhombic tensor
logDiff
DiffFrame
[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
GaussianPP
[GaussianPP LorentzianPP]
lw
GaussianFWHM
[GaussianFWHM LorentzianFWHM]
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-K (λLM,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.
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 handleIf 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
.
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
[evenLmax oddLmax Mmax Kmax]
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
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
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
Opt.pImax
.
GridSize
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'
'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
'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
'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
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'
'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
chili
prints to the command window. If Opt.Verbosity=0
, is is completely silent. 1 prints details about the progress of the computation.
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);
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