Simulation of slow-motion CW-EPR spectra from time-domain trajectories

[B,spc] = cardamom(Sys,Exp)
[B,spc,FID,t] = cardamom(Sys,Exp)
[...] = cardamom(Sys,Exp,Par)
[...] = cardamom(Sys,Exp,Par,Opt)
[...] = cardamom(Sys,Exp,Par,Opt,MD)

cardamom simulates a CW-EPR spectrum using trajectories of orientations, which are obtained from either stochastic dynamics simulations or molecular dynamics (MD) simulations. If stochastic dynamics is to be used, then cardamom uses stochtraj_diffusion or stochtraj_jump to simulate the orientational trajectories required to calculate a cw EPR spectrum. If MD is to be used, then cardamom requires that the information and processed output of the MD simulation (which can itself be done by using mdload) be provided, including the orientational trajectories.

cardamom accepts up to five input arguments

• Sys: static and dynamic properties of the spin system
• Exp: EPR experimental parameters
• Par: simulation parameters for quantum propagation
• Opt: other simulation options
• MD: molecular dynamics parameters

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

Up to four output arguments are returned:

• B: magnetic field axis, in mT
• spc: derivative of the CW-EPR absorption spectrum, arbitrary units
• TDSignal: the time-domain signal, free induction decay (FID)
• t: time axis, in s

The first two output arguments B and spc are mandatory, whereas the last output argument pair TDSignal and t is optional.

Sys is a structure containing the dynamic properties for generating the stochastic rotational trajectories. Many of its fields overlap with those of chili. This allows the user to easily simulate EPR spectra in both frequency and time domains for easy comparison.

Diffusion model

A Brownian diffusion model of rotational motion can be used for local dynamics by declaring Par.Model='diffusion'. For the timescale(s) of motion, one of the fields tcorr, logtcorr, Diff or logDiff should be given. If more than one of these is given, the first in the list logtcorr, tcorr, logDiff, Diff takes precedence over the other(s).

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

It is also possible to specify an orientational potential to simulate restricted diffusion. cardamom accomplishes this by representing the orientational potential function using either a series of Wigner D-functions (see reference for details), a numerical 3D array, or a function handle. To take advantage of one of these feature, use Sys.Potential and assign to it one of the following:

• 2D array (Wigner function expansion). An array of size (N,4), where each row corresponds to a term in the expansion; the first three columns contain the indices L, M, and K; and the fourth column contains the potential coefficients λ^L_M,K (the symbols c and ε are used in the literature as well).

  The ranges of allowed index values are as follows: L≥0; 0≤K≤L; if K=0, then 0≤M≤L; and if K≠0, then -L≤M≤L. Given the associated λ^L_M,K, stochtraj_diffusion will automatically include appropriate λ^L_-M,-K to assure the potential is real-valued. This also means that the input λ^L_0,0 have to be real.

  The full form of the orientational potential is expressed as

  U(Ω) = - k_B T Σ_L,M,Kλ^L_M,K D^L_M,K (Ω),

  where D^L_M,K are the Wigner D-functions, λ^L_M,K are the potential coefficients. The temperature is not needed, as it is divided out in the Langevin equation (which uses only U(Ω)/k_BT).

  For example, to specify an orientational potential using lambda200 to represent the coefficient λ²_0,0 of the commonly used second-order Legendre polynomial, use

  Sys.Potential = [2, 0, 0, 1.2];  % L=2, M=K=0, lambda = 1.2
  

  For details about this type of orientational potential, see K.A. Earle & D.E. Budil, Calculating Slow-Motion ESR Spectra of Spin-Labeled Polymers, in: S. Schlick: Advanced ESR Methods in Polymer Research, Wiley, 2006.

• 3D array (numerical representation). A 3D array of size (N,M,N) representing the orienting potential energy function U(α,β,γ). Dimensions 1, 2, and 3 of this array correspond to the Euler angles α, β, and γ, respectively.

  alpha = linspace(0,2*pi);    % grid over alpha
  beta = linspace(0,pi);       % grid over beta
  gamma = linspace(0,2*pi);    % grid over gamma
  Uarray = lambda*(3*cos(beta).^2-1) % potential energy over grid
  Sys.Potential = Uarray;
  

• Function handle (functional representation). A user-defined function representing U(α,β,γ). This function must by vectorized (accept arrays as inputs) and accept three input arguments, corresponding to the Euler angles α, β, and γ, and yield one output.

  lambda200 = 0.86;
  U = @(a,b,c) lambda200*(3*cos(b).^2-1);
  Sys.Potential = U;
  

Jump model

A jump model with a discrete number of orientational states can be used for local dynamics by declaring Par.Model='jump'. To dictate the inter- and intra-state dynamics, either TransRates or TransProb can be given, with the former taking precedence if both are given.

MD model

An MD model can also be employed, where MD trajectories are used to model local dynamics and obtain additional trajectories. There are three types of MD model that make use of MD trajectories in different ways: Par.Model='MD-direct', where orientational trajectories are directly extracted from the MD trajectory atomic coordinate output; Par.Model='MD-HBD', where the orientations and dynamics of the paramagnetic spin system are used to build and simulate trajectories from a coarse-grained, hindered Brownian dynamics (HBD) model; and Par.Model='MD-HMM', where the dihedral angles of the spin label are projected onto a hidden Markov model (HMM), from which jump trajectories of time-averaged, state-dependent Hamiltonians are simulated.

Par is a structure containing the parameters for generating the stochastic orientational trajectories. The most important field is the model choice:

Three additional mandatory parameters are the spatial dynamics time step dtSpatial, the spin dynamics time step dtSpin, and the number of steps nSteps to use for the calculation of the free-induction decay (FID):

Lastly, the trajectory settings are used to specify the starting orientations for the trajectories using triplets of Euler angles in OriStart = [alpha;beta;gamma], which correspond to the Euler angles α, β, and γ; how many trajectories to simulate, nTraj; and a set of orientations in the lab frame for a powder-like summation, Orients.

If only one set of starting orientations OriStart is given and nTraj is greater than one, then this starting orientation will be repeated for each trajectory. If nTraj is not provided, only one trajectory will be simulated. If OriStart is not provided, then for the case of unrestricted diffusion, a number of starting orientations equal to nTraj will be chosen from a uniform random distribution over the unit sphere, where α, β, and γ are sampled from [0,2π), [0,π), and [0,2π), respectively. If one is simulating using restricted diffusion (Sys.Potential is present), then OriStart will be sampled from the equilibrium distribution corresponding to the potential.

For the lab orientations Orients, if these are not specified but nOrients is declared, then nOrients orientations will be chosen using a spiral grid on the unit sphere. If neither are provided, then the number of lab orientations to sum over will be set equal to Par.nTraj and they will be chosen from a spiral grid.

Opt is a structure with additional simulation options.

MD is a structure with properties and options regarding the usage of MD simulation trajectories to simulate an EPR spectrum.

To simulate a cw EPR spectrum of an ¹⁴N nitroxide spin label with a rotational correlation time of 5 ns using a , use

clear
Sys.g = [2.009, 2.006, 2.002];
Sys.Nucs = '14N';
Sys.A = unitconvert([6,6,36]/10,'mT->MHz');    % hyperfine tensor, G -> MHz
Sys.tcorr = 5e-9;               % 5 ns rotational correlation time

Exp.mwFreq = 9.4;               % microwave frequency, in GHz

Par.Model = 'diffusion';        % rotational diffusion model
Par.dtSpatial = 1e-9;           % spatial propagation time step, in s
Par.dtSpin = 1e-9;              % spin propagation time step, in s
Par.nSteps = 250;               % number of spin propagation time steps

[B,spc] = cardamom(Sys,Exp,Par);

To simulate a CW EPR spectrum using very slow hindered Brownian diffusion, with an orientational potential with coefficient λ²_0,0 = 1.0 and a 250-ns long FID, use

clear
Sys.g = [2.009, 2.006, 2.002];
Sys.Nucs = '14N';
Sys.A = unitconvert([6,6,36]/10,'mT->MHz');    % hyperfine tensor, G -> MHz
Sys.tcorr = 50e-9;  % rotational correlation time, seconds
Sys.Potential = [2, 0, 0, 1.0];

Exp.mwFreq = 9.4; % GHz

Par.Model = 'diffusion';
Par.dtSpatial = 1e-9; % seconds
Par.dtSpin = 1e-9; % seconds
Par.nSteps = 250;
Par.nTraj = 400;

[B,spc] = cardamom(Sys,Exp,Par);

To simulate a cw EPR spectrum using an MD trajectory with a direct model (we assume that the MD trajectory has already been processed and output to the variable MD using mdload), enter

clear
Sys.g = [2.009, 2.006, 2.002];
Sys.Nucs = '14N';
Sys.A = unitconvert([6,6,36]/10,'mT->MHz');    % hyperfine tensor, G -> MHz

Exp.mwFreq = 9.4; % GHz

Par.Model = 'MD-direct';
Par.dtSpin = 1e-9; % seconds, Par.dtSpatial is not needed for MD-direct
Par.nSteps = 250;
Par.nOrients = 400;

[B,spc] = cardamom(Sys,Exp,Par,[],MD);

chili, stochtraj_diffusion, stochtraj_jump, mdload