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. A somewhat easier alternative is to use the graphical user interface that automates many of the steps below. Note that the GUI is several steps behind Spinach kernel development and might not have all the latest features and options exposed.
Contents
Isotopes
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'.
Labels
Optionally, specify a labels 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. 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 zeros 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 zeros 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 zeros 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. 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 shuld be taken to not supply electron coordinates if hyperfine coupling tensors are already specified. |
Example:
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]};
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. |
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 followng 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.
Liquid crystal order matrix
For partially oriented systems, the 3x3 Saupe order matrix may be supplied to enable the simulation of residual dipolar couplings, residual quadrupolar couplings and orientational residuals of all other anisotropic interactions in your system - everything you have specified will be processed. Order matrix specification example with typical values for dilute liquid crystals at room temperature:
inter.order_matrix=diag([1e-3 2e-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 inported 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 parse 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.
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.
Protein and nucleic acid data 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 the J-couplings and some backbone CSA tensors. Particulars and examples are given in the documentation page for protein.m function.
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. Note alsso 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.
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');
Tolerances and algorithm switches
Do not specify any of the parameters below unless you know what you are doing.
If necessary, alter the tolerance parameters for the internal algorithms (sys.tols.*). The defaults are carefully chosen to provide accurate results in every possible case. Performance may often be improved by relaxing some tolerances. Examples:
sys.tols.path_drop=1e-7;
sys.tols.prop_chop=1e-16;
Simulations can vary greatly in their tolerance requirements. Several functions (notably, propagator calculation and optimal control modules) run internal checks and would warn you if they detect any accuracy issues.
If necessary, set the disable switches on the major internal algorithms. The sys.disable variable is a cell array of strings, for example:
sys.disable={'krylov','pt'};
It may be set to disable zero track elimination ('zte'), graph-theoretical analysis ('pt'), permutation symmetry factorization ('symmetry'), Krylov propagation ('krylov'), sparse matrix clean-up ('clean-up') and destination state screening ('dss'). Do not disable any of these algorithms unless you know what you are doing.
If necessary, set the enable switches on experimental algorithms. The sys.enable variable is a cell array of strings, for example:
sys.enable={'gpu','caching'};
It may be set to enable GPU arithmetic ('gpu') and caching of expensive matrix functions ('caching'). On suitable hardware (solid state disks, NVidia Tesla cards), these options can significantly improve performance. Simulations that are constrained by matrix exponentiation stand to benefit particularly. This switch may also be used to force the use of Krylov propagation in the evolution.m module ('krylov').
Revision 3142, authors: Ilya Kuprov, Luke Edwards, Hannah Hogben
