# The spin system

The central concept for spectral simulations with EasySpin is the spin system. It collects all information about the spins in the sample and how they are coupled. In EasySpin, a spin system is represented by a MATLAB structure (the spin system structure) that contains a series of fields holding this information.

The spin system structure contains all information about the spin system and its spin Hamiltonian. Its fields specify spin quantum numbers, interaction parameters, matrices and tensors, relative orientation angles for these matrices and tensors as well as details about broadenings.

The spin Hamiltonian is set up in frequency units (MHz, not angular frequencies), so all the energy parameters of the Hamiltonian have to be given in MHz as well. Hence, for example, the field `A` represents in fact the diagonal of the hyperfine coupling tensor divided by the Planck constant, A/h, and not of A. Remember that all angles are in radians, not in degrees.

Spin system structures are used by many EasySpin functions, among them the simulation functions chili (slow-motion EPR), garlic (isotropic EPR), pepper (solid-state EPR), salt (ENDOR) saffron (pulse EPR), and curry (magnetometry).

A few examples

Before we give a full list of possible fields, here are a couple of examples of spin system definitions. There are two equivalent ways.

```Sys.S = 1/2;
Sys.g = [1.9 2.0 2.1];
Sys.lw = 0.7;            % mT
```
```Sys = struct('S',1/2,'g',[1.9 2.0 2.1],'lw',0.7);
```

Nuclei can be added to a spin system either using a set of fields (`Nucs`, `A`) or using the function nucspinadd.

```Sys.S = 1/2;
Sys.g = [2 2.1 2.2];
Sys.Nucs = '63Cu';
Sys.A = [50 50 460];  % MHz
```
It is more convenient to use nucspinadd:
```Sys = struct('S',1/2,'g',[2 2.1 2.2]);
```

A high-spin Mn2+ system with zero-field splitting is

```Sys = struct('S',5/2,'g',2,'Nucs','55Mn');
Sys.A = -240;   % MHz
Sys.D = 120;    % MHz
```
Spin system parameters

The following groups of parameters can be specified in the spin system structure:

All interaction matrices and tensors can have arbitrary orientations with respect to a molecule-fixed frame of reference (the so-called molecular frame).

Below we list and describe all possible spin system structure fields containing spin Hamiltonian parameters.

The spins

The two fields `S` and `Nucs` are used to specify the electron and nuclear spins constituting the spin system. Both are optional. If `S` is not given, `S=1/2` is assumed. If `Nucs` is not specified, `Nucs=''` is used.

`S`
Gives the electron spin quantum number(s). An arbitrary number of electron spins can be specified. Examples:
```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 `S` is not given, it is automatically set to 1/2.

`Nucs`
A string containing a comma-separated list of nuclear isotopes specifying the nuclear spins in the spin system. An arbitrary number of nuclei can be specified.
```Sys.Nucs = '1H';              % one hydrogen
Sys.Nucs = '63Cu';            % a 63Cu nucleus
Sys.Nucs = '59Co,14N,14N';    % a 59Co and two 14N nuclei
```

If there are no nuclear spins in the system, don't specify this field or set it to `''` (an empty string). Nuclei can be added and removed one at a time using nucspinadd and nucspinrmv.

If not a single isotope, but a natural-abundance mixture of isotopes is needed, just omit the mass number. You can also freely mix single isotopes and natural-abundance mixtures in one spin system.

```Sys.Nucs = 'Cu';          % natural-abundance mixture of 69% 63Cu and 31% 65Cu
Sys.Nucs = 'Cu,14N';      % same as above, plus a pure 14N
Sys.Nucs = 'F,C';         % 100% 19F, plus a mixture of 99% 12C and 1% 13C
```

EasySpin will automatically simulate the spectra of all possible isotopologues (isotopes combinations) and combine them with the appropriate weights. Hyperfine constants and quadrupole couplings are automatically converted between isotopes. The values given by you in the spin system are taken to apply for the most abundant isotope with an appropriate spin (e.g. 63Cu in the case of Cu, 1H for H, 14N of N, 119Sn for Sn).

It is also possible to give custom/enriched isotope mixtures by explicitly giving a list of the mass numbers of all isotopes in parentheses in `Nucs`. Additionally, the abundances should be given in `Abund`. In the simplest case of one atom, this would be

```Sys.Nucs = '(12,13)C';     % for a mixture of 12C and 13C
Sys.Abund = [0.3 0.7];     % 30% of 12C, the rest 13C

Sys.Nucs = '(1,2)H';       % for a mixture of protons and deuterons
Sys.Abund = [0.05 0.95];   % 95% deuterium
```

In the case of multiple atoms, `Abund` should be a cell array with a list of abundances for each atom. For example

```Sys.Nucs = '63Cu,(32,33)S';    % 63Cu with enriched 33S
Sys.Abund = {1,[0.1,0.9]];     % 100% 63Cu, 10% 32S and 90% 33S
```

Custom and natural abundance mixtures and single isotopes can be all used at the same time. Any entry for a natural-abundance atom or for a single isotope in `Abund` is ignored.

```Sys.Nucs = 'Cu,(32,33)S,1H'; % natural Cu with enriched 33S and a pure 1H
Sys.Abund = {1,[0.1,0.9],1}; % one way
Sys.Abund = {[],[0.1,0.9],[]}; % another way, yieldig the same result
```
`Abund`
Specify abundances for custom isotope mixtures. See `Nucs` above.
`n`
Vector of number of equivalent nuclei, e.g.
```Sys.Nucs = '1H,13C';  % one class of 1H and one class of 13C
Sys.n = [2 3];        % 2 protons and 3 carbon-13 spins
```
`Sys.n` can be omitted if all nuclei in `Sys.Nucs` occur only once.

Important: `pepper` does currently not support `Sys.n`. If you have multiple equivalent nuclei, you have to enter each separately.

The g values

The g matrices/tensors for all electron spins in the system are supplied in two fields:

• `g` contains either the 3 principal values of the g tensors, or alternatively the full g tensors.
• `gFrame` contains the Euler angles specifying the orientations of the g tensors relative to the molecular frame.

`g`
Depending on whether the g tensors are rhombic, axial or isotropic, there are different ways you can input the principal values:
• Rhombic: Each row contains the three principal values of the g tensor of one electron spin.

```Sys.g = [2 2.05 2.3];                % rhombic g, for one electron spin
Sys.g = [2 2.1 2.3; 1.9 1.95 2.01];  % rhombic g, for two electron spins
```
• Axial: For axial g tensors, only two values are needed. The first value is the one perpendicular to the unique axis, and the second is the one parallel to it.

```Sys.g = [2.25 2.03];                % axial g, for one electron spin
% = [2.25 2.25 2.03]
Sys.g = [2.25 2.03; 1.99 1.98];     % axial g, for two electron spins
% = [2.25 2.25 2.03; 1.99 1.99 1.98]
```
• Isotropic: If g is isotropic, it is sufficient to give its isotropic value once.

```Sys.g = 2.005;                      % isotropic g, for one electron spin
% = [2.005 2.005 2.005]
Sys.g = [2.0023; 2.0025];           % isotropic g, for two electron spins
% = [2.0023 2.0023 2.0023; 2.0025 2.0025 2.0025]
```
• No value: If no values are given, EasySpin assumes isotropic g values of 2.002319... (see gfree) for all electron spins.

When the principal values of g are given in `g`, the orientation of the g tensor can be specified by Euler angles in `gFrame`, see below.

As an alternative, it is also possible to give the full 3x3 g tensor in `g`. For example

```% full g tensor for one electron spin
Sys.g = [ 2.003397   -0.000431   -0.000004;...
-0.000416    2.003032   -0.000006;...
-0.000014    0.000013    2.002237];

% full g tensors for two electron spins
Sys.g = [2 0 0; 0 2 0; 0 0 2; 2.1 0 0; 0 2.1 0, 0 0 2.1];
```

If you give full g tensors, the Euler angles in `gFrame` are ignored.

`gFrame`
Each row of this array specifies the orientation of the g tensor frame in the molecular frame (see the page on frames). Each row contains the three Euler angles (in radians) for the rotation that transforms the molecular frame to the g tensor frame. If not specified, `gFrame` is assumed to be all-zero, that is, the g tensors of all electron spins are aligned with the molecular frame.
```Sys.gFrame = [0 10 0]*pi/180;                % one electron spin
Sys.gFrame = [0 10 0; 23 -45 67]*pi/180;     % two electron spins
Sys.gFrame = [0 0 0; 0 pi/4 0; 0 -pi/4 0];   % three electron spins
```

With these angles, EasySpin can transform a g tensor from its diagonal eigenframe representation to the molecular frame representation. Here is an explicit example how this is done:

```gpv = [2.1 1.97 2.04];       % principal values
gFrame = [10 34 -2]*pi/180;  % Euler angles, in radians

R_M2g = erot(gFrame)         % transformation matrix
g_g = diag(gpv)              % g in the g frame
g_M = R_M2g.'*g_g*R_M2g      % g in molecular frame
```

which gives the following output

```R_M2g =
0.8220    0.1095   -0.5589
-0.1450    0.9892   -0.0195
0.5507    0.0971    0.8290
g_g =
2.1000         0         0
0    1.9700         0
0         0    2.0400
g_M =
2.0791    0.0154   -0.0278
0.0154    1.9722   -0.0023
-0.0278   -0.0023    2.0587
```

The rows of the transformation matrix `R_M2g` correspond to the g tensor principal axes in molecular frame coordinates. The columns, on the other hand, correspond to the molecular axes in g frame coordinates. See the page on frames for more details.

If a 3x3 transformation matrix `Rg` is given, e.g. from literature, then the corresponding Euler angles can be obtained by

```gFrame = eulang(Rg);
```
`gpa`
This was the way of specifying Euler angles prior to EasySpin 5 and has been superseded by `gFrame`. To convert old `gpa` values to `gFrame`, invert the order and flip the sign of the three Euler angles.
```gpa =    [34  78 -12]*pi/180;    % old form, obsolete
gFrame = [12 -78 -34]*pi/180;    % new form
```

Hyperfine couplings

For each electron-nucleus pair, a hyperfine coupling tensor can be specified. The following fields are used.

`A`
Principal values of the hyperfine interaction matrices A/h, in MHz. Each row contains the principal values of the A tensors of one nucleus. `A(k,:)` refers to nuclear spin `k`. The orientations of the A matrices are given in `AFrame`.
```Sys.A = [-6 12 23];            % 1 electron and 1 nucleus
Sys.A = [10 10 -20; 30 40 50]; % 1 electron and 2 nuclei
```

For axial and isotropic hyperfine tensors, the notation can be shortened, just as in the case of the g tensor.

```Sys.A = [4 10];       % = [4 4 10]  (axial, 1 electron and 1 nucleus)
Sys.A = 34;           % = [34 34 34] (isotropic, 1 electron and 1 nucleus)
Sys.A = [4 10; 1 2];  % = [4 4 10; 1 1 2]  (axial, 1 electron and 2 nucleui)
Sys.A = [7; 3];        % = [7 7 7; 3 3 3] (isotropic, 1 electron and 2 nuclei)
```

If the system contains more than one electron spin, each row contains the principal values of the hyperfine couplings to all electron spins, listed one after the other.

```Sys.A = [10 10 -20 30 40 50];                % 2 electrons and 1 nucleus
Sys.A = [10 10 -20 30 40 50; 1 1 -2 3 4 5];  % 2 electrons and 2 nuclei
```

It is possible to specify full A matrices in `A`. The 3x3 matrices have to be combined like the 1x3 vectors used when only principal values are defined: For different electrons, put the 3x3 matrices side by side, for different nuclei, on top of each other. If full matrices are given in `A`, `AFrame` is ignored.

```Sys.A = [5 0 0; 0 5 0; 0 0 5]                          % 1 electron and 1 nucleus
Sys.A = [[5 0 0; 0 5 0; 0 0 5]; [10 0 0; 0 10 0; 0 0 10]]     % 1 electron and 2 nuclei
Sys.A = [[5 0 0; 0 5 0; 0 0 5], [10 0 0; 0 10 0; 0 0 10]]     % 2 electrons and 1 nucleus
```
`A_`
Spherical form of the hyperfine matrix, in terms of its isotropic, axial and rhombic components. The units are MHz. `A` and `A_` cannot be used at the same time.
```Sys.A_ = 2;         % isotropic component only
Sys.A_ = [2 3]      % isotropic and axial component
Sys.A_ = [2 3 0.4]  % isotropic, axial and rhombic component
```
The cartesian form (as used in `A`) and the spherical form (as used in `A_`) are related by
```% spherical -> cartesian
A(1) = A_(1)-A_(2)-A_(3);
A(2) = A_(1)-A_(2)+A_(3);
A(3) = A_(1)+2*A_(2)

% cartesian -> spherical
A_(1) = mean(A);
A_(2) = (2*A(3)-A(1)-A(2))/6;
A_(3) = (-A(1)+A(2))/2;
```
or if you prefer more compact notation
```A = A_*[1 1 1; -1 -1 2; -1 +1 0];   % spherical -> cartesian
A_ = A*[2 -1 -3; 2 -1 3; 2 2 0]/6;  % cartesian -> spherical
```
For more than one nucleus and more than one electron spin, the array structure is analogous to `A`.
`AFrame`
Array of Euler angles giving the orientations of the various A matrices in the molecular frame, as described above for `gFrame`. If `AFrame` is not specified, it is assumed to be all-zero, that is, all tensors are aligned with the molecular frame. See also erot. Each row of `AFrame` contains the three Euler angles for one nucleus.
```Sys.AFrame = [0 20 0]*pi/180;  % one electron spin, one nucleus
Sys.AFrame = [0 0 0; 0 10 90; 12 -30 34]*pi/180 % one electron spin, three nuclei
```

If there are two or more electron spins, each nucleus has two or more hyperfine tensors, and consequently each row should contain two or more sets of Euler angles.

```Sys.AFrame = [0 20 0, 13 -30 80]*pi/180;               % 2 electrons, 1 nucleus
Sys.AFrame = [0 20 0, 0 0 0; 0 10 0, 0 30 50]*pi/180;  % 2 electrons, 2 nuclei
```

Here is how principal values and angles can be converted to a full matrix:

```Apv = [4 9 13];                   % principal values, MHz
AFrame = [10 45 -30]*pi/180;      % Euler angles
A_A = diag(Apv);                  % full matrix in A frame
R = erot(AFrame);                 % transformation matrix
A_M = R.'*A_A*R;                  % full matrix in molecular frame
```

See the page on frames for more details.

For each nucleus spin, a nuclear quadrupole coupling tensor Q can be given, using the fields `Q` (for the parameters e2qQ/h and η, the principal values, or the full matrix) and `QFrame` (for the orientation).

`Q`
Array containing the quadrupole tensors Q/h for all nuclei, in MHz. There are several possible ways to specify the tensors.
• For axial Q tensors: Specify e2Qq/h for each nucleus
```Sys.Q = 0.7            % one nucleus, eeQq/h = 0.7 MHz
Sys.Q = [0.7, 1.2]     % same for two nuclei
```
• For rhombic Q tensors: Specify e2Qq/h and η for each nucleus, one row with two numbers.
```Sys.Q = [1.2 0.29]         % one nucleus, eeQq/h = 1.2 MHz and eta = 0.29
Sys.Q = [1.2 0.29; 0.1 0]  % same with a second nucleus (second row)
```
• For all Q tensors: specify three principal values of the Q tensor, again one row per nucleus.
```Sys.Q = [-1 -1 2]                % one nucleus, Qxx = -1, Qyy = -1, Qzz = +2 MHz
Sys.Q = [-1 -1 2; 0.2 0.3 -0.5]  % same for two nuclei
```
• You can also give the full 3x3 Q tensor for each nucleus.
```Sys.Q = [-1 0 0; 0 -1 0; 0 0 2]*0.1;   % for one nucleus, diagonal
Sys.Q = [0.125 -0.6495 -1.299; -0.6495 -0.625 0.75; -1.299 0.75 0.5]; % general, one nucleus
Q1 = [-0.9 0 0; 0 -1.1 0; 0 0 2]*0.15;
Q2 = [-0.7 0 0; 0 -1.3 0; 0 0 2]*0.31;
Sys.Q = [Q1; Q2]; % for two nuclei
```
Here is an example of how to convert from e2Qq/h and η to the three principal Q values:
```eeQqh = 1;     % MHz
eta = 0.2;     % unitless
I = 1;         % nuclear spin must be known!
Q = eeQqh/(4*I*(2*I-1)) * [-1+eta, -1-eta, 2]
```

`QFrame`
Euler angles (in radians) describing the orientations of the Q tensors in the molecular frame, analogous to `gFrame` and `AFrame`. `QFrame` should include one row of three angles for each nucleus. The angles are in units of radians, not degrees.

```Sys.QFrame = [0 pi/4 0];                    % one nucleus
Sys.QFrame = [30 45 0; 10 -30 0]*pi/180;    % two nuclei
```

Here is how principal values and angles can be converted to a full matrix:

```Qpv = [-0.9 -1.1 2]*0.125;        % principal values, MHz
QFrame = [10 45 -30]*pi/180;      % Euler angles
R = erot(QFrame);                 % transformation matrix
Q_Q = diag(Qpv);                  % full matrix in Q frame
Q_M = R.'*Q_Q*R;                  % full matrix in molecular frame
```

See the page on frames for more details.

Zero-field splittings

For each electron spin, a zero-field splitting can be specified in the fields `D` and `DFrame`. See also the reference page on the zero-field interaction.

`D`
`D` gives the zero-field splitting tensors for the electron spins in the spin system. It should be in units of MHz (1 cm-1 = 29979 MHz). It can be specified in several different ways.

• Axial: If the zero-field splitting tensor is axial, give the D value for each electron spin. The orientation of the tensor can be specified in the field `DFrame` (see below).

```Sys.D = [200];            % 1 electron spin, D = 200 MHz
Sys.D = [200 -340];       % 2 electron spins, 	D value each, in MHz
Sys.D = [200 -340 1100];  % 3 electron spins, D value each, in MHz
```
• Rhombic: For a zero-field splitting tensor with rhombic asymmetry, give both the D and the E value for each electron spin. The orientation of the tensor can be specified in the field `DFrame` (see below).

```Sys.D = [200 10];            % 1 electron spin, D = 200 MHz and E = 10 MHz
Sys.D = [200 10; 340 90];    % 2 electron spins, D and E value for each spin
```
• Principal values: For both axial and rhombic symmetry, you can give the three principal values of the D tensor. The orientation of the tensor can be specified in the field `DFrame` (see below).

```Sys.D = [-100 -100 200];               % 1 electron spin, Dx = Dy = -100 MHz, Dz = 200 MHz
Sys.D = [-150 -50 200; 450 350 -800];  % 2 electron spins
```
Here is how the principal values can be computed from D and E
```D = 120; E = 15;                     % D and E parameters, in MHz
Sys.D = [-1,-1,2]/3*D + [1,-1,0]*E     % conversion to D principal values
```

It is possible to provide a non-traceless D tensor.

• Full matrix: You can also specify the full D matrix for each electron spin. E.g. for one electron spin

```Sys.D =  [-33.8  -24.1 -122.1; -24.1 -91.2 44.4; -122.1 44.4 125];
```

It is possible to provide a non-traceless D tensor.

Include zeros for any electron spin with S = 1/2.

`DFrame`
nx3 array of Euler angles describing the orientation of the D tensors, completely analogous to `gFrame` and `QFrame`. If absent, it is assumed to be all zeros.

Internally, EasySpin uses the following procedure to compute the full D tensor in the molecular frame from the given principal values and the Euler angles

```Dvals = [-25 -55 80];        % principal values
DFrame = [10 20 0]*pi/180;   % tilt angles
D_D = diag(Dvals);           % full D tensor in its eigenframe

R = erot(DFrame);            % rotation matrix
D_M = R.'*diag(Dvals)*R;     % full D tensor in molecular frame
```

See the page on frames for more details.

High-order zero-field splittings

EasySpin supports a series of high-order electron spin operators in the spin Hamiltonian.

`aF`
Values, in MHz, of the fourth-order parameters a and F used in cubic systems of S=5/2 ions like Fe(III) and Mn(II). For their definition, see the reference page on high-order operators.
```Sys.aF = [10 -3];   % in MHz
```

The operators of both a and F are defined in terms the molecular reference frame, which is assumed to coincide the four-fold symmetry axes of the cubic system.

This can be changed by setting `Sys.aFFrame`. If `Sys.aFFrame=3`, the operator corresponding to a is calculated assuming the reference system coincides with a three-fold symmetry axis of the cubic system. By default, `Sys.aFFrame=4`. As an alternative to `aF`, e.g. for different orientations of the reference frame, the more general high-order parameters `B4` can be used (see below).

`B0, B2, B4, B6, B8, B10, B12`
The coefficients B(k,q) for the extended Stevens operators of the electron spin (see also the function stev). All coefficients should be real, and are understood in units of MHz. The number in the field name indicates the order k of the operator. E.g. B2 is for k = 2, and B4 is k = 4. If any of these fields is not given, it is treated as zero.

Each field that is given should contain a row vector of 2k+1 parameters, in order of decreasing q, running from +k to -k. For example:

```Sys.B2 = [0 0 560 0 0];          % B(2,q) with q = +2,+1,0,-1,-2
Sys.B4 = [0 132 0 0 0 0 0 0 0];  % B(4,q) with q = +4,+3,+2,+1,0,-1,-2,-3,-4
```
The only nonzero elements in the above examples are therefore B(2,0) and B(4,+3).

In the common case that only the q=0 element is needed for a given k, it can be given alone in an abbreviated syntax valid for any k.

```Sys.B2 = [0 0 99 0 0];  % B(2,0) only, full form
Sys.B2 = 99;            % equivalent abbreviated form

Sys.B4 = [0 0 0 0 -8700 0 0 0 0];  % B(4,0) only, full form
Sys.B4 = -8700;                    % equivalent abbreviated form
```

If more than one electron spin is present, specify one row of 2k+1 elements for each electron spin. The first row is for the first electron spin, etc.

```Sys.S = [5/2 2];     % two spins
B20a = 100;          % for the first spin
B20b = -98;          % for the second spin
Sys.B2 = [0 0 B20a 0 0; 0 0 B20b 0 0]; % two rows: one row per spin
Sys.B2 = [B20a; B20b];                 % abbreviated form for two spins: two rows as well
```

Currently, it is not possible to include tilt angles for the principal frames of these high-order interaction tensors. All of them are assumed to be collinear with the molecular frame. By changing the molecular frame, i.e. by tilting all other tensors (g, A, etc), this limitation can be partially circumvented. Still, all the high-order interactions are collinear. Alternatively, you can use `wignerd` to compute a Wigner rotation matrix that can be used to tilt the tensors explicitly.

```angles = rand(1,3)*pi;   % Euler tilt angles, in radians
B2 = [3 4 5 0 2];        % B(2,q) tensor components
R = wignerd(2,angles);   % rotation matrix for rank-2 tensor
B2 = R*B2;               % rotated tensor
Sys.B2 = B2.';
```
Electron-electron spin-spin couplings

For each pair of electron spins, a bilinear coupling matrix (composed of isotropic, anisotropic, and antisymmetric terms) can be given. You can enter it in one of two ways:

• Give the isotropic, antisymmetric and symmetric parts separately in `Sys.J`, `Sys.dvec`, `Sys.eeD`.
• Give the total coupling matrix in `Sys.ee` (with `Sys.eeFrame`).

In addition, it is also possible to specify an isotropic biquadratic exchange coupling for each electron spin pair in `Sys.ee2`. This is not very common.

`J`
1xN array of real

List of isotropic exchange coupling constants (in MHz), one for each pair of electron spins. The associated spin Hamiltonian is +J S1S2 (not -2J S1S2). Here is an example for a two-spin system:

```Sys.S = [3/2 3/2];  % two electron spins
Sys.J = J12;        % coupling constant, in MHz
```

For more than two spins, the pairs are ordered lexicographically, for example for a 4-spin system 1-2, 1-3, 1-4, 2-3, 2-4, 3-4.

```Sys.S = [1/2 1/2 1/2 1/2];          % four spins
Sys.J = [J12 J13 J14 J23 J24 J34];  % six coupling constants, in MHz
```

To convert from the more traditional unit cm-1 to MHz, use

```J_MHz = J_cm*30e3;            % cm^-1 to MHz conversion, approximate
J_MHz = J_cm*100*clight/1e6;  % cm^-1 to MHz conversion, exact
```

You can use either `Sys.J` (together with `Sys.dvec` and `Sys.eeD` if needed) or `Sys.ee` (with `Sys.eeFrame`), but not both at the same time.

`dvec`
Nx3 array of real

List of vectors (one per row, in units of MHz) that describe the antisymmetric part of the coupling between two electron spins (Dzyaloshinskii-Moriya interaction). The associated spin Hamiltonian is d12.(S1xS2), where d12 is the interaction vector.

```Sys.S = [1/2 1/2];    % two spins
d12 = [dx dy dz];     % units: MHz
Sys.dvec = d12;       % one interaction vector
```

For more than two spins, the vectors are ordered in lexicographic order, for example 1-2, 1-3, 2-3 for a three-spin system.

```Sys.S = [1 1 1];        % three spins
d12 = [0 0 1.4e4];      % units: MHz
d23 = [0 0 2.1e4];
d13 = [0 0 3.7e4];
Sys.J = [d12;d13;d23];  % three interaction vectors, one per row
```

You can use either `Sys.J` (together with `Sys.dvec` and `Sys.eeD` if needed) or `Sys.ee` (with `Sys.eeFrame`), but not both at the same time.

`eeD`
Nx3 array of real

List of vectors (one per row, in units of MHz) that describe the symmetric part of the coupling between two electron spins (dipolar interaction). Each row contains the three principal values of the symmetric interaction tensor. Their sum has to give zero (i.e. the tensor has to be traceless).

```Sys.S = [1 3/2];           % two spins
Sys.eeD = [-1 -1 2]*140;   % principal value of dipolar coupling tensor, MHz
```

For more than two spins, the vectors are ordered in lexicographic order, for example 1-2, 1-3, 2-3 for a three-spin system.

```Sys.S = [3/2 3/2 3/2];    % three spins
D12 = [-1 -1 2]*1.4e4;    % units: MHz
D23 = [-1 -1 2]*2.1e4;
D13 = [-1 -1 2]*3.7e4;
Sys.J = [D12; D13; D23];  % three tensors
```

You can use either `Sys.J` (together with `Sys.dvec` and `Sys.eeD` if needed) or `Sys.ee` (with `Sys.eeFrame`), but not both at the same time.

`ee`
Nx3 or 3Nx3 array of real

Principal value of the electron-electron interaction matrices, in MHz. Each row corresponds to the diagonal of an interaction matrix (in its eigenframe). For two electron spins, `ee` contains one row only.

```Sys.S = [1/2 1/2];     % two electron spins
Sys.ee = [50 50 100];  % principal values of one coupling matrix, MHz
```

For more than two electron spins, the pairs are lexicographically ordered according to the indices of the electron spins involved. If n is the number of electron spins, there are N = n(n-1)/2 rows. For example, for 4 spins there are 6 rows with the principal values for the interaction of spins 1-2, 1-3, 1-4, 2-3, 2-4, 3-4, in this order.

```Sys.S = [1/2 1/2 1/2];        % three spins
ee12 = [50 50 100];
ee13 = [10 10 -30];
ee14 = [10 10 -30];
ee23 = [0 0 0];
ee24 = [0 0 0];
ee34 = [80 80 80];
Sys.ee = [ee12; ee13; ee14; ee23; ee24; ee34];  % 6 coupling matrices, MHz
```

If only isotropic couplings are needed, use `Sys.J` (see above).

It is also possible to specify the full 3x3 interaction matrices instead of the 3 principal values and 3 Euler angles. These matrices combine the isotropic, antisymmetric and symmetric parts of the interaction. For a 2-electron system, `ee` is a single 3x3 array.

```Sys.S = [1/2 1/2];                   % two spins
Sys.ee = [50 0 0;0 50 0; 0 0 100];   % one coupling matrix, MHz
```

For more electrons, the 3x3 matrices are stacked on top of each other, to give a 3Nx3 array.

```Sys.S = [1/2 1/2 1/2];                   % three spins
ee12 = [50 0 0; 0 50 0; 0 0 100];
ee13 = diag([10 10 -30]);
ee23 = zeros(3);
Sys.ee = [ee12; ee13; ee23];             % three coupling matrices
```

Note that you can use either `Sys.ee` (together with `Sys.eeFrame`) or `Sys.J` (together with `Sys.dvec` and `Sys.eeD` if needed), but not both at the same time.

`eeFrame`
Nx3 array containing the Euler angles describing the orientation of the electron-electron interaction matrices (specified in `ee`) in the molecular frame. Each row contains the Euler angles for the corresponding row in `ee`.

This is fully analogous to `AFrame` (see above). Check out the page on frames for more details on frames.

`ee2`
Contains the isotropic biquadratic exchange coupling, in MHz, for each pair of electron spins. The spin pairs are ordered as in `ee`.
```Sys.S = [3/2 3/2];          % two electron spins
Sys.ee2 = 130;              % one biquadratic coupling 1-2

Sys.S = [3/2 3/2 3/2];      % three electron spins
Sys.ee2 = [130 150 190];    % couplings 1-2, 1-3, 2-3

Sys.S = [1 1 1 1];               % four electron spins
Sys.ee2 = [130 0 130 130 0 130]; % couplings 1-2, 1-3, 1-4, 2-3, 2-4, 3-4
```

EasySpin defines the biquadratic exchange term in the spin Hamiltonian as +j(S1.S2)2, with a plus sign. In the literature, it is sometimes defined with a negative sign, so be careful when using literature values.

Nuclear-nuclear spin-spin couplings

For each pair of nuclear spins, a bilinear coupling matrix (composed of isotropic, anisotropic, and antisymmetric terms) can be given. This term is specified in the field `Sys.nn`, with the orientation in `Sys.nnFrame`.

In practice, this term is essentially always too small to be of any relevance for EPR or ENDOR spectra. EasySpin provides it so its effects on spectra can be explored.

`nn`
Nx3 or 3Nx3 array of real

Principal values of the nuclear-nuclear interaction matrices, or full interaction matrices, in MHz. For providing principal values, each row contains the three principal values of one interaction matrix. For two nuclear spins, `nn` contains one row only.

```Sys.Nucs = '1H,1H';        % two nuclear spins
Sys.nn = [10 20 30]*1e-6;  % principal values of coupling matrix between spins 1 and 2, MHz
```

For more than two nuclear spins, the pairs are lexicographically ordered according to the indices of the nuclear spins involved. If n is the number of nuclear spins, there are N = n(n-1)/2 rows. For example, for 4 spins there are 6 rows with the principal values for the interaction of spins 1-2, 1-3, 1-4, 2-3, 2-4, 3-4, in this order.

```Sys.Nucs = '1H,1H,19F';        % three spins
nn12 = [50 50 100]*1e-6;
nn13 = [10 10 -30]*1e-6;
nn14 = [10 10 -30]*1e-6;
nn23 = [0 0 0];
nn24 = [0 0 0];
nn34 = [80 80 80]*1e-6;
Sys.nn = [nn12; nn13; nn14; nn23; nn24; nn34];  % 6 coupling matrices, MHz
```

If the principal values are given, the orientation of the interaction matrices can be specified via Euler angles in the field `nnFrame` (see below).

It is also possible to specify the full 3x3 interaction matrices instead of the 3 principal values and 3 Euler angles. For two nuclei, `nn` is a single 3x3 array.

```Sys.Nucs = '1H,1H';                       % two nuclear spins
Sys.nn = [50 0 0;0 50 0; 0 0 100]*1e-6;   % one coupling matrix, MHz
```

For more nuclei, the 3x3 matrices are stacked on top of each other, to give a 3Nx3 array.

```Sys.S = '1H,1H,1H';                   % three spins
nn12 = [50 0 0; 0 50 0; 0 0 100]*1e-6;
nn13 = diag([10 10 -30])*1e-6;
nn23 = zeros(3);
Sys.nn = [nn12; nn13; nn23];          % three coupling matrices, MHz
```
`nnFrame`
Nx3 array of real

Array containing the Euler angles describing the orientation of the nuclear-nuclear interaction matrices (specified in `nn`) in the molecular frame. Each row contains the three Euler angles for the corresponding row in `nn`.

This is fully analogous to `eeFrame` (see above). Check out the page on frames for more details on frames.

Spin-orbit interaction
For each electron spin an orbital angular momentum can be defined in the field `L`. It is described by the crystal field splitting, the orbital reduction factor and the spin-orbit coupling. Currently, spin-orbit interaction is supported by the simulation functions pepper (solid-state EPR)and curry.

Please cite Joscha and friends Journal of unsolved questions, 2020, 87, 92-109 when using spin-orbit interaction.

`L`
Defines the orbital angular momenta to which the electron spins should be coupled. `L` is an optional field. If given it need to have the same number of entries as `S` and the entries have to be positive integer or 0. Examples:
```Sys = struct('S',[1/2 1/2], 'L',[2,2]); % two spin-1/2, each associated to L = 2
Sys = struct('S',[1/2 3/2], 'L',[0,1]); % one spin-1/2 without L and a spin 3/2 with L = 1
Sys = struct('S',1, 'L',3);             % one spin 1 associated to a L = 3
```
`CF0, CF2, CF4, CF6, CF8, CF10, CF12`
The crystal field of the orbital angular momentum L is in EasySpin described by the coefficients CF(k,q) for the extended Stevens operators (see also the function stev). They are analog to the coefficients B(k,q) for the electron spin (see Higher-order operators), but acting on the orbital angular momentum. All coefficients should be real, and are understood in units of MHz. The number in the field name indicates the order k of the operator. E.g. `CF2` is for k = 2, and `CF4` is k = 4. If any of these fields is not given, it is treated as zero. `L` must be given, otherwise the CF are ignored.

Each field that is given should contain a row vector of 2k+1 parameters, in order of decreasing q, running from +k to -k. For example:

```Sys.CF2 = [0 0 560 0 0];          % CF(2,q) with q = +2,+1,0,-1,-2
Sys.CF4 = [0 132 0 0 0 0 0 0 0];  % CF(4,q) with q = +4,+3,+2,+1,0,-1,-2,-3,-4
```
The only nonzero elements in the above examples are therefore CF(2,0) and CF(4,+3).

In the common case that only the q=0 element is needed for a given k, it can be given alone in an abbreviated syntax valid for any k.

```Sys.CF2 = [0 0 99 0 0];  % CF(2,0) only, full form
Sys.CF2 = 99;            % equivalent abbreviated form

Sys.CF4 = [0 0 0 0 -8700 0 0 0 0];  % CF(4,0) only, full form
Sys.CF4 = -8700;                    % equivalent abbreviated form
```

If more than one electron spin (and therefore more than one orbital angular moment) is present, specify one row of 2k+1 elements for each orbital angular moment. The first row is for the first orbital angular moment, etc.

```Sys.S = [5/2 2];     % two spins
Sys.L = [1 2];       % with two orbital angular momenta
CF20a = 100;          % for the first orbital angular moment
CF20b = -98;          % for the second orbital angular moment
Sys.CF2 = [0 0 CF20a 0 0; 0 0 CF20b 0 0]; % two rows: one row per orbital angular moment
Sys.CF2 = [CF20a; CF20b];                 % abbreviated form for two orbital angular moment, two rows as well
```

Currently, it is not possible to include tilt angles for the principal frames of these crystal field interaction tensors. All of them are assumed to be collinear with the molecular frame. By changing the molecular frame, i.e. by tilting all other tensors (g, A, etc), this limitation can be partially circumvented. Still, all the crystal-field interactions are collinear. Alternatively, you can use `wignerd` to compute a Wigner rotation matrix that can be used to tilt the tensors explicitly.

```angles = rand(1,3)*pi;   % Euler tilt angles, in radians
CF2 = [3 4 5 0 2];        % CF(2,q) tensor components
R = wignerd(2,angles);   % rotation matrix for rank-2 tensor
CF2 = R*CF2;               % rotated tensor
Sys.CF2 = CF2.';
```
`orf`
Orbital reduction factor entering in the Zeeman Hamiltonian and the spin-orbit interaction for each orbital angular momentum. `orf` is an optional field, if omitted it is assumed as 1 for all orbital angular momenta. The number of elements in `orf` has to match the number of elements in `L`.
```Sys.orf = 3/2;          % orf of a single effective l =1 as calculated by Lines for octahedral Co(II)
Sys.orf = [0.97, 0.93]; % orf for two orbital angular momenta
```
`soc`
Spin-orbit coupling in MHz. It can also contain higher orders of the spin-orbit coupling, as sometimes used for heavier elements. Each electron spin can only be coupled to one orbital angular momentum. If more than one electron is present, the first row contain the coupling between the first electron spin and the first orbital angular momentum and so on. Examples:
```% three electron spins, with linear soc parameter of 170, 250 cm-1 and 0:
Sys.soc = [170; 250; 0]*clight*1e-4;
% as above, plus the first spin has a cubic soc parameter of 10 cm-1:
Sys.soc = [170 0, 10; 250, 0, 0; 0, 0, 0]*clight*1e-4;

Sys = struct('S',[1/2 1], 'L',[2,1]); % two electron spin with S =1/2 and 1 and L = 2 and 1
Sys.soc = [200; 100]*clight*1e-4;
% linear soc parameter between S = 1/2 and L = 2 is 200 cm-1
% linear soc parameter between S = 1 and L = 1 is 100 cm-1
```
General Parameters for Spin Hamiltonian expanded in Magnetic Field and Electron Spin

The Spin-Hamiltonian is in general a polynom of the magnetic field and the spin operator. This general expansion is implemented for electron spins (not for nuclei). Currently, general parameter are supported by the simulation functions pepper (solid-state EPR)and curry. The coefficients can be provided in the following way:

`Ham`
The coefficients C(pB, pS, p, q) should be real and in units of MHz (mT)-pB.The pB, pS, p are integrated in the field name in this order. If any of these fields are not given, they are treated as zero. Each field that is given should contain a row vector of 2p+1 parameters, in order of decreasing q, running from +p to -p. For example:
```Sys.Ham134 = [0 0 36 0 0 0 12 0 0 0]; % pB = 1, pS = 3, p = 4, q = +4,+3,+2,+1,0,-1,-2,-3,-4
```

```%convert the more common D and E:
Sys.Ham022 = [sqrt(2)*E 0 sqrt(2/3)*D 0 0];
%convert an isotropic g value:
Sys.Ham110 = -bmagn/planck*1e-9*sqrt(3)*g;
```
There is a number of fields by which line broadenings can be specified. `lw` and `lwEndor` are line widths (FWHM) which are used for convolution of the simulated spectrum. All others are so-called strains and describe Gaussian distributions in the associated spin Hamiltonian parameters.