Spin system specification

Any Spinach calculation must begin with a call to create.m and basis.m constructor functions. They process spin system and simulation formalism specifications and create the spin_system object – the primary data structure that is used to store spin system information in Spinach. The first function call is:

    spin_system=create(sys,inter);

Follow the instructions below to prepare sys and inter arguments.

Isotopes and text labels

Specify the spin system composition by giving a list of isotope names, for example:

    sys.isotopes={'1H','1H','19F','235U'};

An electron may be selected by specifying 'E'. Electron shells with multiplicity higher than 2 may be selected by specifying a multiplicity after the electron label, e.g. 'E3' requests a spin-1 electron. All known isotopes are supported, you can specify your own particles by editing the spin.m file in the kernel. Spinach supports particles of any spin, including spin zero. Spin zero particles only have the unit state – they do not influence matrix dimension and all couplings to them are effectively ignored. To request a spin zero particle, simply specify the corresponding isotope, for example '12C'. Spinach also has a spin-zero "ghost" particle built in. A ghost particle does not participate in any dynamics and may be specified using 'G'.

Optionally, specify a label for each spin by giving a list of strings, for example:

    sys.labels={'CA','CB','HB2','HB3'};

Labels are printed next to spin interaction summaries – this makes them easier to read for large spin systems. When they are unique, labels may be used (via idxof.m function) to refer to spins by their labels. Labels are also used by protein NMR spectroscopy modules to identify different types of atoms – when a dedicated protein pulse sequence (such as hncoca.m) is run, these labels must be set to the standard PDB atom identifiers. If this field is omitted, all labels are set to empty strings. PDB and BMRB import function protein.m sets these labels automatically.

Zeeman interactions

The primary magnet field must be specified in units of Tesla. Example:

    sys.magnet=14.1;

For simulations of systems that are not inside any permanent magnet (such as those encountered in RYDMR spectroscopy) you must specify zero.

If your system has chemical shifts or g-tensors present, the corresponding interactions for all particles in the system may be specified as scalars, 3x3 matrices, or eigenvalues + Euler angles (in radians). If multiple specifications are supplied, they are added together.

Zeeman interaction specification
Variable name	Variable type	Content
inter.zeeman.eigs, inter.zeeman.euler	[1 × nspins] cell arrays of [1 × 3] row vectors	Eigenvalues of chemical shift tensors (in ppm for nuclei) and g-tensors (in Bohr magneton units for electrons) with Euler angles (in radians). Individual cells in the array may be left empty, in which case pure magnet frequencies determined by the free particle magnetogyric ratio are assumed.
inter.zeeman.matrix	[1 × nspins] cell array of [3 × 3] matrices	Full chemical shift tensors (in ppm for nuclei) and g-tensors (in Bohr magneton units for electrons) as matrices. Individual cells in the array may be left empty, in which case pure magnet frequencies determined by the free particle magnetogyric ratio are assumed.
inter.zeeman.scalar	[1 × nspins] cell array of real numbers	Isotropic chemical shifts (in ppm for nuclei) and g-factors (in Bohr magneton units, for electrons). Individual cells in the array may be left empty, in which case pure magnet frequencies determined by the free particle magnetogyric ratio are assumed.

Examples:

    inter.zeeman.eigs={[7 15 -22] ...
                       [11 18 -29]};
    inter.zeeman.euler={[pi/5 pi/3 pi/11] ...
                        [pi/6 pi/7 pi/15]};

    inter.zeeman.matrix={[5 0 0; 0 5 0; 0 0 5] ...
                         [5 0 0; 0 5 0; 0 0 5] ...
                         [2.0023 0 0; 0 2.0025 0; 0 0 2.0027]};

    inter.zeeman.scalar={2.0023 1.0 2.0 3.0};

Spin-spin interactions

If you have spin-spin couplings in your system, they may be specified as scalars, 3x3 matrices, or eigenvalues + Euler angles, or (for giant spin Hamiltonian parameters) irreducible spherical tensor coefficients. If multiple specifications are supplied, they are added together.

Spin-spin coupling specification
Variable name	Variable type	Content
inter.coupling.eigs, inter.coupling.euler	[nspins × nspins] cell array of [1 × 3] matrices	Eigenvalues of coupling tensors (in Hz) with Euler angles (in radians). Bilinear coupling is introduced by specifying a coupling between two different spins. Quadratic coupling (e.g. quadrupolar) is introduced by specifying a coupling between a spin and itself. Individual cells in the array may be left empty, in which case zeros are assumed.
inter.coupling.matrix	[nspins × nspins] cell array of [3 × 3] matrices	Full coupling tensors as matrices (in Hz). Each element of the cell array is account-ed for, so the couplings must be divided by two if a symmetric cell array is supplied. Individual cells in the array may be left empty, in which case zeros are assumed.
inter.coupling.scalar	[nspins × nspins] cell array of reals	Isotropic couplings (in Hz). Individual cells in the array may be left empty, in which case zeros are assumed.
inter.coordinates	[nspins × 1] cell array of [1 × 3] row vectors	Cartesian coordinates of every spin (in Angstroms), used to determine point dipolar interactions. If a cell corresponding to a particular spin is left empty, that spin is assumed to not have any dipolar interactions with the rest of the system. Care should be taken to not supply nuclear coordinates if hyperfine coupling tensors are already specified for those nuclei.
inter.giant.coeff	A [1 x nspins] cell array of cell arrays with a three-element vector followed by a five-element vector, followed by a seven-element vector and so on	Coefficients (in the order of descending projection quantum number) in front of the three first-rank irreducible spherical tensors, five second-rank irreducible spherical tensors, seven third-rank irreducible spherical tensors, and so on - in the giant spin Hamiltonian model. Stevens operator coefficients must be translated into irreducible spherical tensor operator coefficients before use. Spinach does not natively support Stevens operators, and never will.
inter.giant.euler	[1 × nspins] cell array of cell arrays of [1 × 3] row vectors	Euler angles (in radians) specifying the molecular frame orientation of the interactions specified via the irreducible spherical tensor coefficients given in inter.giant.coeff field.

Examples:

    inter.coupling.scalar={0 50; 0 0};
    inter.coupling.eigs{2,2}=[1e4 1e4 -2e4];
    inter.coupling.euler{2,2}=[0 0 0];
    inter.coordinates={[0.0 0.0 0.0] ...
                       [0.0 0.0 1.5]};

    inter.giant.coeff={{[0 0 0],[0 0 -4.65e8 0 0],[1e7 0 0 2e7 0 0 1e7]}};
    inter.giant.euler={{[0 0 0],[0 0 0],[0 0 0]}};

Spin-spin interactions may be specified in a variety of equivalent ways. The table below provides suggestions on specifying all common interactions.

Ways of specifying spin interactions
Nuclear chemical shift	Use inter.zeeman.scalar for isotropic chemical shifts, inter.zeeman.matrix for anisotropic chemical shift tensors supplied as matrices, or inter.zeeman.eigs & inter.zeeman.euler for anisotropic chemical shift tensors specified as eigenvalues and Euler angles. Mind the difference between chemical shielding and chemical shift - Spinach requires chemical shift.
Inter-nuclear J-coupling	Use inter.coupling.scalar; couplings that are specified multiple times, e.g. between spin 1 and 2, and then again between spin 2 and 1, will be added together.
Inter-nuclear dipolar coupling	Use inter.coordinates if nuclear coordinates are known (they will be converted into a dipolar interaction matrix internally), or inter.coupling.matrix for dipolar coupling supplied as a matrix (make sure the matrix is traceless), or inter.coupling.eigs & inter.coupling.euler for dipolar interactions supplied as eigenvalues and Euler angles (make sure the eigenvalues sum up to zero). When supplying dipolar interactions in any way other than Cartesian coordinates, make sure that the matrices you have supplied actually correspond to a feasible arrangement of particles in 3D space. You would likely want to set the cut-off distance beyond which the dipolar couplings should be ignored to simplify and accelerate the simulation. This is done by specifying the corresponding tolerance. The default is 100 Angstrom.
Nuclear quadrupolar coupling	Best specified as an interaction of the nucleus with itself. Use inter.coupling.matrix (make sure the matrix is traceless), or inter.coupling.eigs & inter.coupling.euler for quadrupolar interactions specified as eigenvalues and Euler angles (make sure the eigenvalues sum up to zero). Nuclear quadrupolar interaction conventions in the literature are sometimes obscure, eeqq2nqi.m function might help.
Electron g-factor or g-tensor	Use inter.zeeman.scalar for isotropic g-factors, or inter.zeeman.matrix for anisotropic g-tensors supplied as matrices, or inter.zeeman.eigs & inter.zeeman.euler for anisotropic g-tensors specified as eigenvalues and Euler angles. Be careful with the numerical accuracy of the values that you provide; at least six decimal places are recommended.
Electron-nuclear hyperfine coupling	Use inter.coupling.scalar for isotropic hyperfine interactions, or inter.coupling.matrix for hyperfine tensors specified as matrices, or inter.coupling.eigs & inter.coupling.euler for hyperfine tensors specified as eigenvalues and Euler angles. If the hyperfine coupling is purely dipolar and occurs at a distance of more than 10 Angstrom, electron and nuclear locations may be specified using inter.coordinates, they will be converted into a hyperfine coupling automatically. In this case you would likely want to set the cut-off distance beyond which the dipolar couplings should be ignored to simplify and accelerate the simulation. This is done by specifying the corresponding tolerance. The default is 100 Angstrom.
Electron zero-field splitting	Best specified as an interaction of the electron with itself. Use inter.coupling.matrix (make sure the matrix is traceless), or inter.coupling.eigs & inter.coupling.euler for ZFS specified as eigenvalues and Euler angles (make sure eigenvalues sum up to zero).
Inter-electron exchange coupling	Use inter.coupling.scalar; couplings that are specified multiple times, e.g. between spin 1 and 2, and then again between spin 2 and 1, will be added together.
Inter-electron dipolar coupling	Use inter.coordinates if the point dipole approximation is good enough and electron coordinates are known (they will be converted into a dipolar interaction matrix internally), or inter.coupling.matrix for dipolar couplings supplied matrices (make sure the matrices are traceless), or inter.coupling.eigs & inter.coupling.euler for dipolar interactions supplied as eigenvalues and Euler angles (make sure the eigenvalues sum up to zero). When supplying dipolar interactions in any way other than Cartesian coordinates, make sure that the matrices you have supplied actually correspond to a feasible arrangement of particles in 3D space. You would likely want to set the cut-off distance beyond which the dipolar couplings should be ignored to simplify and accelerate the simulation. This is done by specifying the corresponding tolerance. The default is 100 Angstrom.
Giant spin Hamiltonian terms	Use inter.giant.coeff and inter.giant.euler; remember that coefficients refer to the irreducible spherical tensors (in the order of descending projection quantum number), not Stevens operators. Stevens operator coefficients must be translated before use.

Periodic boundary conditions

Periodic boundary conditions in one, two, or three dimensions may be specified by supplying one, two, or three lattice translation vectors, in Angstroms, as shown in the following examples:

    inter.pbc={[0 0 10]};

    or

    inter.pbc={[0 0 10],[5 5 0]};

    or

    inter.pbc={[0 0 10],[5 5 0],[10 0 0]};

Lattice translation vectors need not be orthogonal, but must be linearly independent. Periodic boundary conditions are used for the calculation of dipolar couplings. The number of images beyond which the couplings are to be ignored may be specified by setting the corresponding tolerance. The default value is to account for two images in each direction. Note also that the default proximity cut-off beyond which dipolar couplings are ignored is 100 Angstroms. This may also be changed by altering the tolerance settings.

Magnetic susceptibility tensors

For systems with a rapidly relaxing paramagnetic centres, the 3x3 magnetic susceptiblity tensors (in cubic Angstroms) may be provided as a cell array of 3x3 matrices, for example:

    inter.suscept.chi={[0.0883   -0.0904    0.0822
                       -0.0904   -0.1011   -0.0149
                        0.0822   -0.0149    0.0128]};

Locations of the susceptibility tensors must be provided as a cell array of Cartesian coordinates in Angstroms, for example:

    inter.suscept.xyz={[1.0 2.5 3.9]};

Specifying susceptibility tensors triggers the calculation of paramagnetic shift tensors on all nuclei. These tensors are added to the CSAs and automatically enter Redfield theory when it is used, meanging that Curie relaxation and its cross-correlations are accounted for automatically. Note that only spherical rank 0 and 2 components are processed at the moment; the antisymmetric contribution to the paramagnetic tensor is currently ignored.

Liquid crystal order matrix

For partially oriented systems, the 3x3 Saupe order matrices (one per chemical component, as a cell array of 3x3 matrices) may be supplied to enable the simulation of residual dipolar couplings, residual quadrupolar couplings, and orientational residuals of any other anisotropic interactions in your system. Order matrix specification example with typical values for dilute liquid crystals at room temperature:

    inter.order_matrix={diag([1e-3 2e-3 -3e-3])};

For a system with two chemical compartments:

    inter.order_matrix={diag([1e-3 2e-3 -3e-3]), diag([2e-3 1e-3 -3e-3])};

At simulation time, the transformation of interaction tensors into their partially averaged from is handled by liquid.m context function. It may be performed manually by calling residual.m function.

Gaussian and ORCA import

Magnetic interaction parameters and atomic coordinates may be imported directly into sys and inter data structures from Gaussian03/09 and ORCA logs. In both cases, the log is first parsed using either gparse.m or oparse.m function and then the data is imported into Spinach by calling g2spinach.m function, for example:

    % Parse a DFT calculation log
    props=gparse('../standard_systems/alanine.log');

    % Import data into Spinach
    [sys,inter]=g2spinach(props,{{'C','13C'},{'N','15N'}},[182.1 264.5],[]);

Further details on the parameters and options for gparse.m, oparse.m and g2spinach.m are given in the manual pages of those functions. Note that chemical shielding and chemical shift are different things.

Some examples of Gaussian inputs: gaussian_a.txt, gaussian_b.txt, gaussian_c.txt - note that the extra printing flag (#p) should always be present.

For chemical shift tensors, the anisotropies computed with electronic structure methods (DFT, MP2, etc.) are usually satisfactory for the purposes of relaxation theory, but the isotropic chemical shifts can differ significantly from the experimental values. After the import is complete, the isotropic parts of chemical shift tensors may be replaced with their experimental values using shift_iso.m function.

To verify the chemical shifts after the import and a call to create.m, either inspect the console output, or use chemshifts.m function. Imports from multiple electronic structure theory calculations may be merged using merge_inp.m function.

SIMPSON import

Spin system information may also be read from the spinsys{} field of a SIMPSON *.in file using the following syntax:

    [sys,inter]=s2spinach(filename);

The function absorbs isotope information, all Zeeman tensors and all coupling tensors, and returns Spinach sys and inter data structures with that information. Structures from multiple imports may be merged using merge_inp.m function.

Protein and nucleic acid import

Protein spin system composition and interaction information may be loaded from a pair of protein database files – a PDB file with atomic coordinates and a BMRB file with chemical shifts. The following call:

    [sys,inter]=protein(pdb_file,bmrb_file,selection)

would automatically form the necessary data structures, estimate all J-couplings and some backbone CSA tensors. Particulars and examples are given in the documentation page for protein.m function. Note that both PDB and BMRB data formats are hard to parse; please refer to the protein NMR example files to see the particular formats that Spinach is able to parse.

Nucleic acid data may be imported in a similar way:

    [sys,inter]=nuclacid(pdb_file,shift_file,options)

The the manual page of nuclacid.m for further information. Spinach example set contains several examples of RNA NMR simulations.

Console and file output control

Spinach prints large amounts of information to Matlab console. All console output with the exception of fatal error messages may be redirected to an ASCII text file by setting

    sys.output='filename';

or suppressed altogether (useful for large sets of simulations and parameter fitting runs) by setting

    sys.output='hush';

The default is to print all output to Matlab console. An exception is powder averages inside powder.m, singlerot.m and doublerot.m, where the console output is suppressed unless the user specifically indicates that it shouldn't be by setting

    parameters.verbose=1;

in the experiment parameters. Note also that slave processes in parallel calculations are not permitted to write their output to the log because the result of multiple processes writing to a file simultaneously is usually a mess. However, all threads are allowed to print to the console, so do not expect the output of highly parallel runs to be easily readable.

Spinach has no way of telling when your calculation is finished and therefore would not close the log file for you. It is therefore the user's responsibility to close all open files at the end of the calculation, for example by calling

    fclose('all');

This command closes all open files.

Tolerances and algorithm switches

The default internal accuracy settings in Spinach are very tight and intended to prevent any unpleasant surprises in all circumstances. In each individual case, however, these tight settings might be unnecessary. The tolerance parameters that you should consider relaxing from their default values are:

    sys.tols.prox_cutoff  - the distance beyond which dipolar couplings are ignored, default is 100 Angstrom

    sys.tols.inter_cutoff - the amplitude below which spin-spin interactions are ignored, default is 1e-10 rad/s

For example, in liquid state protein systems there are, in general, millions of dipolar interactions, but only those over distances shorter than about 5 Angstrom are in practice important. Setting the proximity cut-off tolerance to 5 Angstrom reduces the number of dipolar interactions to less then a hundred thousand. This accelerates the simulatios dramatically.

Another important aspect to consider is that the internal heuristics in Spinach is optimised for very large spin systems. If your system has fewer than five spins, there is no point running sophisticated trajectory-level state space analysis. In that case, some performance may be gained by disabling it:

    sys.disable={'trajlevel'};

Spinach plays it safe with memory requirements and would automatically switch to slower (but much less memory-hungry) Krylov propagation if it sees a large Liouvillian. To force the faster matrix exponential path, specify:

    sys.disable={'krylov'};

and so on - multiple internal algorithms may be disabled in this way, see the expert level options page for further information.

There are also a few internal algorithms that are very useful, but for now considered experimental, such as GPU support and matrix exponential caching. To enable them, set:

    sys.enable={'gpu','caching'};

On suitable hardware (solid state disks, NVidia Tesla cards), these options can significantly improve performance. See the expert level options page for further information.

Version 2.8, authors: Ilya Kuprov, Luke Edwards, Hannah Hogben