stochtraj_diffusion

Simulation of stochastic rotational trajectories.

Syntax

stochtraj_diffusion(Sys)
stochtraj_diffusion(Sys,Par)
stochtraj_diffusion(Sys,Par,Opt)
[t,RTraj] = stochtraj_diffusion(...)
[t,RTraj,qTraj] = stochtraj_diffusion(...)

Description

stochtraj_diffusion computes stochastic trajectories of rotational diffusion using orientational quaternions.

stochtraj_diffusion accepts up to three input arguments

Sys: dynamic properties
Par: (optional) simulation parameters for monte carlo integrator and trajectories
Opt: (optional) simulation options

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

Up to three output arguments are returned:

t: time axis, in seconds
RTraj: trajectory in terms of rotation matrices
qTraj: trajectory in terms of quaternions

If no output argument is given, stochtraj_diffusion plots one of the trajectories.

Input: System dynamics

Sys is a structure containing the dynamic properties for generating the stochastic 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 comparison.

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

tcorr

Rotational correlation time, in seconds.

1 number: isotopic diffusion
2 numbers [txy tz]: anisotropic diffusion with axial diffusion tensor
3 numbers [tx ty tz]: anisotropic diffusion with rhombic diffusion tensor

For example,

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

Instead of tcorr, Diff can be used, see below. If tcorr is given, Diff is ignored. The correlation time tcorr and the diffusion rate Diff are related by tcorr = 1/(6*Diff).

logtcorr

Base-10 logarithm of the correlation time, offering an alternative way to input the correlation time. If given, tcorr, logDiff and Diff are ignored.
Use this instead of tcorr for least-squares fitting with esfit.

Diff

Rotational diffusion rates (principal values of the rotational diffusion tensor), in seconds^-1.

one number: isotopic diffusion tensor
two numbers: input [Dxy Dzz] gives axial tensor [Dxy Dxy Dzz]
three numbers: rhombic tensor [Dxy Dxy Dzz]

Diff is ignored if logtcorr, tcorr or logDiff is given.

logDiff

Base-10 logarithm of Diff. If given, Diff is ignored.

Potential

It is also possible to specify an orienting potential energy function to simulate restricted rotational diffusion. stochtraj_diffusion accomplishes this by representing the orienting 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 orienting 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;
```

Input: Simulation parameters

Par is a structure containing the parameters for calculating the stochastic trajectories. The parameters can be categorized into two groups: Monte Carlo settings and trajectory settings.

The Monte Carlo settings specify the parameters used to perform the time integration in the simulation. One can provide an array of time points t along which to integrate, a number of steps nSteps and a time step dt, or a number of steps nSteps and a total time for simulation tmax.

t

An array of time points along which to integrate, in seconds. If this is provided, no other Monte Carlo settings need to be specified.

nSteps

The total number of time steps used to perform integration.

dt

The time step of the integration, in seconds.

tmax

The total time of the simulation, in seconds.

Lastly, the trajectory settings are used to specify the number of trajectories to be simulated, nTraj, and the starting states of each trajectory using OriStart.

OriStart

Starting orientations for trajectories, given as an array of size (3,1), (1,3), or (3,nTraj) containing triplets of the Euler angles α, β, γ, in radians. If OriStart is not provided, then a number of starting orientations equal to nTraj will be chosen from a uniform random distribution, where α, β, and γ are sampled from [0,2π), [0,π), and [0,2π), respectively.

If simulating restricted diffusion using an orienting potential, it is strongly recommended to not specify Par.OriStart. In this case, it is best to let stochtraj_diffusion determine the optimal distribution of starting orientations for the trajectories to achieve faster convergence, which is dictated by sampling the Boltzmann distribution associated with the orienting potential.

nTraj

The total number of trajectories to simulate.

If only one starting orientation is given and nTraj is greater than one, then this starting orientation will be used for each trajectory. If OriStart is not provided, then starting orientations will be chosen uniformly at random over the unit 4-sphere. If nTraj is not provided, only one trajectory will be simulated.

Input: Simulation options

Opt is a structure with additional simulation options.

checkConvergence

0 (default), 1
If equal to 1, after the first nSteps of the trajectories are calculated, stochtraj_diffusion will check for both inter- and intra-trajectory convergence using the Gelman-Rubin R statistic (see A. Gelman and D. Rubin, Statistical Science 7, 457 (1992) for details). The trajectories are considered converged if R < 1+Opt.convTolerance. If this condition is not satisfied after the nSteps, then propagation will be extended by either a length of time equal to the average correlation time or by 20% of the current trajectory length, whichever is larger. If convergence is not obtained, then propagation will be extended by either a length of time equal to the average of tcorr or by 20% more time steps, whichever is larger.

convTolerance

Convergence tolerance for Gelman-Rubin R statistic. The threshold for R is 1 + Opt.convTolerance, e.g. if Opt.convTolerance = 1e-6, then the threshold R is 1.000001.

Verbosity

0 (default), 1
Determines how much information stochtraj_diffusion prints to the screen. If Opt.Verbosity=0, it is completely silent. 1 prints details about the progress of the computation.

Example

To simulate a 200-ns long, fast-motion Brownian diffusion trajectory starting from the north pole on the unit sphere, enter

Sys.tcorr = 5e-9;       % 5 ns rotational correlation time

Par.dt = 1e-9;          % 1 ns time step
Par.nSteps = 1e3;       % 1000 time steps
Par.OriStart = [0,0,0]; % start at north pole

[t,RTraj] = stochtraj_diffusion(Sys,Par);

To simulate 400 1-μs long, slow-motion trajectories in the presence of an orienting potential with coefficient λ²_0,0 = 1.2, and with starting orientations chosen randomly, use

Sys.tcorr = 50e-9;  % 50 ns rotational correlation time
Sys.Potential = [2, 0, 0, 1.2];  % L=2, M=K=0, lambda = 1.2

Par.dt = 1e-9;       % 1 ns time step 
Par.nSteps = 1000;   % 1000 steps
Par.nTraj = 400;     % 400 trajectories

[t,RTraj] = stochtraj_diffusion(Sys,Par);

Algorithm

stochtraj_diffusion solves the stochastic rotational Langevin equation by numerically propagating quaternions. The quaternions specify orientations of the body-fixed coordinate frame with respect to the space-fixed frame. Since both isotropic and anisotropic diffusion tensors are diagonal in the body-fixed frame, defining the quaternions in this manner considerably simplifies the problem.

For full details of the methods used here, see:

D. Sezer, J. H. Freed, and B. Roux, J. Chem. Phys. B 128, 165106 (2008).