Kernel functions

Spinach kernel functions provide basic simulation infrastructure - operators, states, propagators, and other simulation building blocks. They also perform various operations that may be viewed as elementary - time propagation, coherence selection, rotating frame transformations, etc.

All kernel functions require spin system information to be provided, as a spin_system object in the first argument:

    answer=kernel_function(spin_system,...);

This page is an overview of kernel functions - click on function names to see further particulars. Kernel function code is clear, compact and extensively commented. It is usually a good idea to take a look inside a function before using it.

Simulation assumptions

assume.m

The composition of the spin Hamiltonian depends on the assumptions that are made about the spin system. Spinach requires the user to specify the assumptions (using assume.m function) before hamiltonian.m function can be permitted to run, for example:

    spin_system=assume(spin_system,'deer');

This command sets up interaction specifications inside the spin_system object in such a way that a call to hamiltonian.m would produce a Hamiltonian in the rotating frame that is typically used in DEER spectroscopy. See the description of assume.m for further information on the assumption sets supported by Spinach out of the box. You can create your own assumption sets by editing that function.

Perturbative corrections to rotating frame transformations should be obtained by setting the assumptions to 'labframe' and then calling rotframe.m to perform the rotating frame transformation numerically to the perturbation theory order that you require. More specialised interaction representation transformations (for example, with respect to zero-field splitting) should be applied to the laboratory frame Hamiltonian using intrep.m function.

residual.m

In partially ordered systems, such as liquid crystals, this function triggers the application of partial averaging to the anisotropic interactions. When this function is invoked, all interactions in the spin_system structure are overwritten with their partial averaging residuals.

Order matrix must be supplied at the spin system specification stage, otherwise this function would return an error.

This function must be called after the relaxation superoperator is generated because it destroys interaction anisotropy information.

Elementary operators

hamiltonian.m

Hamiltonian and its rotational decomposition. In Liouville space formalisms (sphten-liouv, zeeman-liouv) a Hamiltonian superoperator (left product, right product or commutation) is returned. In Hilbert space (zeeman-hilb formalism) the Hamiltonian matrix is returned. The matrix produced by the function is always sparse. A typical call for the isotropic part of the Hamiltonian commutation superoperator in Liouville space is:

    spin_system=assume(spin_system,'nmr');
    H=hamiltonian(spin_system,'comm');

A typical call to obtain the full anisotropic Hamiltonian at a specific system orientation is:

    spin_system=assume(spin_system,'nmr');
    [H,Q]=hamiltonian(spin_system,'comm');
    H=H+orientation(Q,[alpha beta gamma]);

This function should never be called repeatedly at every orientation in a powder average - a much more efficient procedure exists whereby hamiltonian.m is called once with two outputs before the powder loop is executed. At each instance of the loop, these outputs are supplied to orientation.m function that returns the Hamiltonian at that specific orientation, for example:

    spin_system=assume(spin_system,'nmr');
    [I,Q]=hamiltonian(spin_system,'comm');
    <begin powder loop>
         H=I+orientation(Q,[alpha beta gamma]);
         <use the Hamiltonian for something>
    <end powder loop>

See the description of hamiltonian.m for further information. Hamiltonian generation is parallelized.

relaxation.m

This function generates relaxation superoperators in Liouville space calculations (sphten-liouv and zeeman-liouv formalism). The call format is:

    R=relaxation(spin_system,[alpha beta gamma]);

The second parameter only applies to the H-strain relaxation model. All relaxation theory parameters are set during the spin system and interaction specification that is supplied to create.m function. The relaxation superoperator is generated as a negative definite matrix that should be added to the Hamiltonian commutation superoperator in the following way:

    L=H+1i*R;

See the spin system specification section and the description of relaxation.m for further details. Relaxation superoperator generation is parallelized.

kinetics.m

This function generates chemical kinetics superoperators in Liouville space calculations (sphten-liouv and zeeman-liouv formalism). The call format is:

    K=kinetics(spin_system);

All chemical kinetics parameters are set during the spin system and interaction specification that is supplied to create.m function. The kinetics superoperator should be added to the Hamiltonian commutation superoperator in the following way:

    L=H+1i*K;

See the spin system specification section and the description of kinetics.m for further details.

propagator.m

This function bulds exponential propagators using scaled and squared power series, trying, insofar as possible, to keep the matrix sparse. The function requires the Liouvillian matrix and the time step:

    P=propagator(spin_system,L,timestep);

If GPU support has been switched on during the call to create.m (see the spin system specification section), the system has a card that supports CUDA and the matrix meets the dimensionality requirements, then the evaluation is carried out on the GPU.

If expensive function caching has been switched on during the call to create.m (see the spin system specification section), and the matrix is big enough to bother, the function would inspect its cache and either store the output if the argument has not been encountered before, or retrieve the result without re-computing it if the argument is found in the cache.

GPU and caching options are particularly recommended for calculations that are dominated by matrix exponentiation step, such as some solid state NMR simulations.

operator.m

This function generates user-specified operators and superoperators. Two types of syntax are supported. To generate a sum of specified operators on all spins of a specified type, specify the spin and the operator as character strings, for example:

    sum_Lz=operator(spin_system,'Lz','13C');

To generate a specific product operator or superoperator, use cell arrays, for example:

    LzSp=operator(spin_system,{'Lz','L+'},{1,2});

A more detailed description is available in operator.m page.

orientation.m

This function uses the isotropic Hamiltonian I and the spherical components of the anisotropic Hamiltonian Q returned by the hamiltonian.m function, for example:

    [I,Q]=hamiltonian(spin_system);

to produce the Hamiltonian at the specified spin system orientation:

    H=I+orientation(Q,[alpha beta gamma]);

A call to orientation.m is much faster than re-creating the Hamiltonian from scratch. It should always be used in powder averages and other operations in which the system Hamiltonian at multiple orientations is required.

Elementary states

state.m

Generates user-specified state vectors (in Liouville space) and density matrices (in Hilbert space). Two types of syntax are supported. To generate a sum of specified states on all spins of a specified type, specify the spin and the state as character strings, for example:

    sum_Lz=operator(spin_system,'Lz','13C');

To generate a specific product state, use cell arrays, for example:

   LzSp=operator(spin_system,{'Lz','L+'},{1,2});

A more detailed description is available in state.m page.

equilibrium.m

Returns thermal equilibrium state vectors (in Liouville space) and density matrices (in Hilbert space). The temperature is specified during the call to create.m, see the spin system specification section. The call syntax in Hilbert space is:

    rho=equilibrium(spin_system,H,Q,euler_angles);

where the last two arguments are optional and only needed when the system is anisotropic. In Lioville space, make sure that the Hamiltonian supplied to the function is a laboratory frame left side product superoperator:

    [H,Q]=hamiltonian(assume(spin_system,'labframe'),'left');

Sending a commutation superoperator or a rotating frame Hamiltonian would produce incorrect results.

singlet.m

Returns two-spin singlet state. The numbers of the spins involved should be supplied to the function:

    rho=singlet(spin_system,spin_a,spin_b)

A state vector is returned in Liouville space and a density matrix in Hilbert space.

triplet.m

Returns the components of the two-spin triplet state. The numbers of the spins involved should be supplied to the function:

    [Tp,T0,Tm]=triplet(spin_system,spin_a,spin_b)

State vectors are returned in Liouville space and density matrices in Hilbert space.

Time evolution

Different functions should be used depending on the nature of the time propagation problem:

evolution.m

Long-range evolution under a static Liouvillian, including detection and refocussing stages of NMR and EPR pulse sequences. This function also handles the time integrals of observables, such as those occurring in spin chemistry calculations.

krylov.m

Same functionality as evolution.m with the difference that Krylov propagation is used without any attempt to reduce the problem dimension. This module should be used with very large (dimension >> 10,000) Liouvillians that cannot be reduced or exponentiated without overflowing the system memory. Krylov method avoids all matrix-matrix operations and only uses sparse matrix-vector products for propagation. This saves significant amounts of memory, but may be slower then running evolution.m.

step.m

This is a subset of krylov.m that is optimized for one-off events such as hard pulses. This function should be used for sliced propagation under time-dependent Liouvillians.

isergen.m

Iserles product quadrature generators. These are product equivalents of trapezium and Simpson quadrature that are much more accurate than the usual piecewise-constant propagator approximation.

iserstep.m

Iserles solvers for the Lie equation, where the generator can depend on the state. These are product equivalents of trapezium and Simpson quadrature that are much more accurate than the usual piecewise-constant propagator approximation.

Shaped pulses and gradients

grad_pulse.m

Analytical approximation to a gradient pulse - works by wiping the populations of states with non-zero Zeeman frequency. To model gradient dynamics faithfully with a possibility of refocussing, use Fokker-Planck formalism within imaging.m context.

grad_sandw.m

Analytical approximation to a gradient sandwich - works by wiping the populations of states with non-zero effective Zeeman frequency. To model gradient dynamics faithfully with a possibility of refocussing, use Fokker-Planck formalism within imaging.m context.

cartesian2polar.m

Translates Cartesian coordinates and derivatives with respect to those coordinates into the polar form.

chirp_pulse_af.m

Chirp pulse waveform with a sine bell amplitude envelope, amplitude-frequency representation.

chirp_pulse_xy.m

Chirp pulse waveform with a sine bell amplitude envelope, Cartesian representation.

pmlg5.m

PMLG5 phase cycle.

polar2cartesian.m

Translates polar coordinates and derivatives with respect to those coordinates into the Cartesian form.

pulse_shape.m

Shaped pulse waveforms.

read_wave.m

Reads JCAMP-DX waveform files.

restrans.m

Resonator response function convolution.

rsequence.m

Symmetry-based R-sequences.

rseq_compiler.m

R-sequence compiler.

save_wave.m

Writes a pulse waveform to a JCAMP-DX pulse file.

sawtooth.m

Returns a sawtooth waveform.

shaped_pulse_af.m

Shaped pulse in amplitude-frequency coordinates using Fokker-Planck formalism.

shaped_pulse_xy.m

Shaped pulse function using Cartesian coordinates.

spinal.m

SPINAL phase cycle.

triwave.m

Returns a triangular waveform.

vg_pulse.m

Veshtort-Griffin shaped pulses.

wave_basis.m

Common basis sets for the expansion of pulse waveforms.

Trajectory and state analysis

There are two functions that return trajectory content and similarity information based on the basis set data available within sphten-liouv formalism.

trajan.m

Performs trajectory analysis by partitioning the population of the system state vector into user-specified subspaces and plotting the populations of those subspaces as a function of time. This function is useful in checking the accuracy of restricted basis sets -- if significant population is detected in the highest available spin correlation order, the basis is likely to be too small.

trajsimil.m

Performs similarity analysis of two spin system trajectories by partitioning the dynamics into user-specified subspaces and plotting the difference in the population of those subspaces as a function of time.

stateinfo.m

Prints the state vector norm and the list of the most populated basis states in the order of decreasing population. This function requires a spherical tensor basis set.

Coherence order selection

At the state vector level, coherence order and correlation order selection are performed analytically by classifying basis states. This procedure is only available in sphten-liouv formalism.

coherence.m

Performs analytical coherence selection by computing spin state coherence orders from the basis descriptor and zeroing the population of undesired states. Example:

    rho=coherence(spin_system,rho,{{'13C',[1 -1]},{'1H',-1}})

This command keeps the population of the states that have coherence order ((1 or -1 on 13C) and (-1 on 1H)) and zeros the population of all other states.

It is also possible to select coherence orders during basis set construction, so that irrelevant states are never considered in the simulation to begin with -- see the correlation and coherence order screening section of the basis set specification page.

correlation.m

Performs analytical correlation order selection by computing spin state correlation orders from the basis descriptor and zeroing the population of undesired states. Example:

    rho=correlation(spin_system,rho,[1 2],'13C')

This command keeps correlation orders 1 and 2 involving 13C nuclei.

It is also possible to select correlation orders during basis set construction, so that irrelevant states are never considered in the simulation to begin with -- see the correlation and coherence order screening section of the basis set specification page.

homospoil.m

Emulates the effect of a homospoil gradient pulse by identifying the states that would be defocussed by the pulse and zeroing their amplitudes. This function uses the basis set information and is therefore only available within sphten-liouv formalism.

Rotating frames and average Hamiltonians

Rotating frames can either be applied analytically at the Hamiltonian construction time (see assume.m) or numerically after the Hamiltonian is built (in that case assume.m should be called with a 'labframe' argument). The following functions are relevant:

carrier.m

The carrier Hamiltonian, containing all Zeeman terms prescribed by the current magnet and magnetogyric ratios:

    H=carrier(spin_system,spins)

average.m

Average Hamiltonian theory over the period of some periodic process. This module should be used when the details of the coupling Hamiltonian (i.e. the positive and the negative frequency parts under the periodic process in question) are known:

    H=average(spin_system,Hp,H0,Hm,omega,theory)

This function supports Waugh, Krylov-Bogolyubov and brute-force period averaging methods.

rotframe.m

Numerical rotating frame transformation module that supports perturbative corrections to the rotating frame Hamiltonian:

    H=rotframe(spin_system,C,H,isotope,order)

This function can be expensive and should only be used when analytical rotating frame transformation is not available for some reason.

frqoffset.m

In high-field NMR simulations where perturbative corrections to the rotating frame are not necessary, this function is used to shift the frequency offset position for a specific set of spins:

    H=frqoffset(spin_system,H,parameters)

where the parameters data structure contains information about spins and their frequency offsets.

Decoupling and spin locking

decouple.m

This function obliterates all interactions and populations in the subspace of states that involve user-specified spins in any way:

    [L,rho]=decouple(spin_system,L,rho,spins)

The specified spins would not contribute to the system dynamics until the Liouvillian is rebuilt from scratch. This is "ideal" decoupling. To model the process explicitly, make your Hamiltonian time-dependent.

spinlock.m

This function obliterates all spin-spin correlations and all magnetization components other than those along the indicated direction:

    rho=spinlock(spin_system,Lx,Ly,rho,direction)

This is an "ideal" spin lock. To model the process explicitly, add RF terms to the Hamiltonian.

Version 2.6, authors: Ilya Kuprov