The spin system

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

Sys.S = 1/2;
Sys.g = [2 2.1 2.2];
Sys.Nucs = '63Cu';
Sys.A = [50 50 460];  % MHz

Sys = struct('S',1/2,'g',[2 2.1 2.2]);
Sys = nucspinadd(Sys,'63Cu',[50 50 460]);

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

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 electron spin
Sys.S = [1/2, 1/2];     % two S=1/2 electron spins
Sys.S = [1, 1, 1/2];    % two S=1 electron spins and one S=1/2 electron spin

If S is not specified, it is automatically set to 1/2.

A character array 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 character array). 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 in the spin system are taken to apply for the most abundant isotope with an appropriate spin (e.g. Cu: 63Cu for both gn and quadrupole moment; H: 1H for g value, 2H for quadrupole moment; N: 14N for both; Sn: 119Sn for both; the only three cases where the reference isotopes for gn and quadrupole moment differ are H, Xe and Hg, for all others either the most abundant isotope has spin 1 or larger [N], or there is no isotope with spin 1 or larger [C]).

It is also possible to specify custom isotope mixtures by listing the mass numbers of all isotopes in parentheses in Nucs and giving the associated abundances 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];   % 5% 1H, 95% 2H

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

Specify abundances for custom isotope mixtures. See Nucs above.

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 not support Sys.n. If you have multiple equivalent nuclei, enter each separately.

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.

Alternatively, it is possible to provide 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, EasySpin uses them as provided, i.e. the tensors are not symmetrized, and the Euler angles in gFrame are ignored.

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 from molecular frame to g frame
R_g2M = R_M2g.';             % transformation matrix from g frame to molecular frame
g_g = diag(gpv)              % g in g frame
g_M = R_g2M*g_g*R_g2M.'      % 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);

Prior to EasySpin 5, the field gpa was used instead of gFrame. To convert old gpa values to gFrame, invert the order and flip the sign of the three Euler angles.

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

In the case where the associated nucleus is a natural or custom isotope mixture, the given hyperfine value refers to the I>=1/2 isotope with the highest abundance. See the description under Sys.Nucs.

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.

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 A matrix in A frame
R_M2A = erot(AFrame);             % transformation matrix from molecular frame to A frame
R_A2M = R_M2A.';                  % transformation matrix from A frame to molecular frame
A_M = R_A2M*A_A*R_A2M;            % full A matrix in molecular frame

See the page on frames for more details.

Prior to EasySpin 5, the field Apa was used instead of AFrame. To convert old Apa values to AFrame, invert the order and flip the sign of the three Euler angles.

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 e²Qq/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 e²Qq/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 e²Qq/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]

And here is a demonstration of how to convert from the three principal values of Q to e²Qq/h and η:

Q = [-0.5 -1.0 1.5]       % MHz
I = 1;                    % nuclear spin must be known!

eeQqh = 2*I*(2*I-1)*Q(3)  % MHz
eta = (Q(1)-Q(2))/Q(3)    % unitless

See also the reference page on nuclear quadrupole coupling.

In the case where the associated nucleus is a natural or custom isotope mixture, the given quadrupole data refers to the I>=1 isotope with the highest abundance. See the description under Sys.Nucs.

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_M2Q = erot(QFrame);             % transformation matrix from molecular frame to Q frame
R_Q2M = R_M2Q.';                  % transformation matrix from Q frame to molecular frame
Q_Q = diag(Qpv);                  % full matrix in Q frame
Q_M = R_Q2M*Q_Q*R_Q2M.';          % full matrix in molecular frame

See the page on frames for more details.

Prior to EasySpin 5, the field Qpa was used instead of QFrame. To convert old Qpa values to QFrame, invert the order and flip the sign of the three Euler angles.

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; 300 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.

D_ gives the zero-field splitting tensors for the electron spins in the spin system defined in terms of D and the ratio E/D. D should be in units of MHz (1 cm^-1 = 29979 MHz), E/D is unitless.

Sys.D_ = [200 0.05];             % 1 electron spin, D = 200 MHz and E/D = 0.05, i.e. E = 10 MHz
Sys.D_ = [200 0.05; 300 0.3];    % 2 electron spins, D and E/D value for each spin

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_M2D = erot(DFrame);              % transformation matrix from molecular frame to D frame
R_D2M = R_M2D.';                   % transformation matrix from D frame to molecular frame
D_M = R_D2M*diag(Dvals)*R_D2M.';   % full D tensor in molecular frame

See the page on frames for more details.

Prior to EasySpin 5, the field Dpa was used instead of DFrame. To convert old Dpa values to DFrame, invert the order and flip the sign of the three Euler angles.

The coefficients B(k,q) for the extended Stevens operators of the electron spin (see also the function stev). All coefficients must be real-valued, and are understood in units of MHz. The number in the field name indicates the rank k of the tensor 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).

Positive (and zero) q correspond to to the cosine tesseral operators, and negative q correspond to the sine tesseral operators.

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.

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

To include tilt angles for Sys.B1, Sys.B2, etc, use the fields Sys.B1Frame, Sys.B2Frame, etc. See below.

Both Sys.B2 and Sys.D provide parameters for the rank-2 zero-field splitting term in the spin Hamiltonian. Therefore, only one of the two can be used. The two spin systems Sys1 and Sys2 in the following example give identical Hamiltonians:

D = 1000; % MHz
E = D*0.2; % MHz
ang = [pi/3 pi/7 pi/11]; % tilt angles
Sys1.D = [D E];
Sys1.DFrame = ang;
Sys2.B2 = [E 0 D/3 0 0]; 
Sys2.B2Frame = ang;

In the special case of high-spin 3d⁵ ions such as Fe(III), Mn(II), and Cr(I) in axially distorted cubic environments, the fourth-order terms are often written in terms of the traditional a and F parameters instead of B₄. For their definition, see the reference page on high-order operators. EasySpin versions prior to 6 had a dedicated field Sys.aF for these. Since EasySpin 6, this is no longer supported. Instead, the fourth-order terms should be specified using Sys.B4 and Sys.B4Frame. Pay attention that the frames for the a and the F term are different.

3-element arrays containing the Euler angles describing the orientation of the tensors specified in B1, B2, etc. (see above) in the molecular frame. Use B1Frame for the rank-1 tensor in B1, B2Frame for the rank-2 tensor in B2, and so forth.

These fields are fully analogous to AFrame (see above). See the page on frames for more details on frames.

Here is an example for rank 2:

Sys.B2 = [0 30 100 0 0];   % in MHz
Sys.B2Frame = [0 pi/3 0];  % in radians

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 S₁S₂ (not -2J S₁S₂). 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 = unitconvert(J_cm,'cm^-1->MHz');            % cm^-1 to MHz conversion

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

Nx1, Nx2, or Nx3 array of real

List of values that describe the symmetric part of the coupling between two electron spins (dipolar interaction), in units of MHz. Each row corresponds to one electron-electron pair. It can contain 1, 2, or 3 numbers.

• 1 number: This gives the axial dipolar coupling.
• 2 numbers: The first is the axial dipolar coupling, and the second one is the rhombic dipolar coupling.
• 3 numbers: They give the three principal values of the symmetric dipolar interaction tensor. The trace of a dipolar tensor is zero.

Here are a few examples:

Sys.S = [1/2 1/2];    % two spins

% The following two inputs give the same axial dipolar tensor
Sys.dip = 23;          % axial dipolar coupling, MHz
Sys.dip = 23*[1 1 -2]; % principal values of axial dipolar coupling tensor, MHz

% The following two inputs give the same rhombic dipolar tensor
Sys.dip = [23 3];                    % axial and rhombic dipolar coupling, MHz
Sys.dip = 23*[1 1 -2] + 3*[1 -1 0];  % principal values of axial dipolar coupling tensor, MHz

The axial and rhombic couplings and the principal values of the dipolar tensor are related via

% axial and rhombic -> principal values
dip_pv = dip_axial*[1 1 -2] + dip_rhombic*[1 -1 0];

% principal values -> axial and rhombic
dip_rhombic = (dip_pv(1)-dip_pv(2))/2
dip_axial = ((dip_pv(1)+dip_pv(2))/2-dip_pv(3))/3

For more than two spins, specify the couplings for electron pairs 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
dip12 = 1.4e4;    % units: MHz
dip23 = 2.1e4;
dip13 = 3.7e4;
Sys.dip = [dip12; dip13; dip23];  % three tensors

To specify the orientation of the dipolar tensor, use Sys.eeFrame.

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

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 d₁₂.(S₁xS₂), where d₁₂ 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.dvec = [d12;d13;d23];  % three interaction vectors, one per row

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

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 = [80 80 20];  % 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 = [80 80 20];
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 = [80 0 0;0 80 0; 0 0 20];   % 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 = [80 0 0; 0 80 0; 0 0 20];
ee13 = diag([10 10 -30]);
ee23 = zeros(3);
Sys.ee = [ee12; ee13; ee23];             % three coupling matrices

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

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.

Prior to EasySpin 5, the field eepa was used instead of eeFrame. To convert old eepa values to eeFrame, invert the order and flip the sign of the three Euler angles.

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)², with a plus sign. In the literature, it is sometimes defined with a negative sign, so be careful when using literature values.

Nx3 or 3Nx3 array of real

Principal values of the CS tensors, or full CS tensors. When providing principal values, each row contains the three principal values of one CS tensor.

Sys.Nucs = '1H';             % one nuclear spin
Sys.sigma = 1-[4 5 9]*1e-6;  % principal values of shielding tensor

Nx3 array of real

Array containing the Euler angles describing the orientation of the chemical shielding tensors (specified in sigma) in the molecular frame. Each row contains the three Euler angles for the corresponding row in sigma.

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

Nx3 or 3Nx3 array of real

Principal values of the nucleus-nucleus 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

Nx3 array of real

Array containing the Euler angles describing the orientation of the nucleus-nucleus 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.

Gives the quantum numbers for orbital angular momenta to which the electron spins are coupled. L is an optional field. If given, it must have the same number of entries as S and the entries have to be nonnegative integers. The kth orbital angular momentum is coupled to the kth spin angular momentum. Examples:

% one spin 1 associated to a L = 3
Sys.S = 1;
Sys.L = 3;

% two spin-1/2, each associated to L = 2
Sys.S = [1/2 1/2];
Sys.L = [2 2];

% one spin-1/2 without L and a spin 3/2 with L = 1
Sys.S = [1/2 3/2];
Sys.L = [0 1];

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

Orbital g values entering in the Zeeman Hamiltonian, one for each orbital angular momentum in L. gL is an optional field, if omitted it is assumed equal to 1 for all orbital angular momenta. The number of elements in gL must match the number of elements in L.

Sys.gL = 3/2;          % gL for a single orbital angular momentum
Sys.gL = [0.97 0.93];  % gL for two orbital angular momenta

Spin-orbit coupling constants, 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; 

% two electron spin with S =1/2 and 1 and L = 2 and 1
Sys.S = [1/2 1];
Sys.L = [2 1];
Sys.soc = [200; 100]*clight*1e-4; 
% linear soc of 200 cm-1 between S = 1/2 and L = 2 
% linear soc of 100 cm-1 parameter between S = 1 and L = 1

The coefficients C(p_B, p_S, p, q) should be real and in units of MHz (mT)^-p_B.The p_B, p_S, 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;

Sys.initState = {state, basis}

  En(1) = (1/3)*D - E; % Tx
  En(2) = (1/3)*D + E; % Ty
  En(3) = -(2/3)*D;    % Tz