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

Rotational correlation time, in seconds.

• 1 number txyz: isotopic diffusion, corresponding to [txyz txyz txyz]
• 2 numbers [txy tz]: anisotropic diffusion with axial diffusion tensor, corresponding to [txy txy tz]
• 3 numbers [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

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

Principal values of the rotational diffusion tensor, in s^-1 (equivalently, rad² s^-1).

• one number Dxyz: isotopic diffusion tensor, corresponding to [Dxyz Dxyz Dxyz]
• two numbers [Dxy Dzz]: axial tensor, corresponding to [Dxy Dxy Dzz]
• three numbers [Dxy Dxy Dzz]: rhombic tensor

Base-10 logarithm of Diff. Use this instead of Diff for least-squares fitting with esfit.

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.

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).

• 1 element: GaussianPP
• 2 elements: [GaussianPP LorentzianPP]

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).

• 1 element: GaussianFWHM
• 2 elements: [GaussianFWHM LorentzianFWHM]

An array of coefficients for the orientational potential, with four numbers per row: L, M, K, and λ. The potential is expressed as U(Ω) = - k_B T Σ_L,M,Kλ^L_M,KD^L_M,K, where the D^L_M,K are Wigner functions and λ^L_M,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 (λ^L_M,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:

• If the environment of cylindrically symmetric (e.g. in a uniaxial liquid crystal), then use a potential with M=0 terms.
• If the tumbling paramagnetic molecules is cylindrically symmetric, use a potential with K=0 terms only.

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.

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.

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 cos²β - 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.

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.

-1 (default) or +1
Determines which electron spin transitions to include in the spin basis. If set to false, includes all transitions with p_S = m_S'-m_S'' = -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 p_S = +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.

vector, one element per nucleus
Determines for each nucleus the maximum nuclear coherence order |Δm_I| to include in the spin basis. 0 corresponds to only allowed EPR transitions, 1 includes additionally forbidden EPR transitions with |Δm_I|=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.

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.

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.

'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.

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.

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.

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.

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.

'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.

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.

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);