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.

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 equivlent ways. The table below provides suggestions on specifying all common interactions.

Ways of specifying magnetic interactions
Nuclear chemical shift	inter.zeeman.scalar for isotropic chemical shifts, inter.zeeman.matrix for anisotropic chemical shift tensors supplied as matrices, 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	inter.coupling.scalar
Nuclear dipolar coupling	inter.coordinates if nuclear coordinates are known (they will be converted into a dipolar interaction matrix internally), inter.coupling.matrix for dipolar coupling supplied as a matrix (make sure it's traceless), inter.coupling.eigs & inter.coupling.euler for dipolar interactions supplied as eigenvalues and Euler angles (make sure eigenvalues sum up to zero).
Nuclear quadrupolar coupling	inter.coupling.matrix (as an interaction of the nucleus with itself) or inter.coupling.eigs & inter.coupling.euler for quadrupolar interactions specified as eigenvalues and euler angles (make sure 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	inter.zeeman.scalar for isotropic g-tensors, inter.zeeman.matrix for anisotropic g-tensors supplied as matrices, inter.zeeman.eigs & inter.zeeman.euler for anisotropic g-tensors specified as eigenvalues and Euler angles.
Hyperfine coupling	inter.coupling.scalar for isotropic hyperfine couplings, inter.coupling.matrix for hyperfine tensors specified as matrices, 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 coordinates may be specified using inter.coordinates, they will be converted into a hyperfine coupling automatically. Note that sys.tols.prox_cutoff tolerance may need to be increased from its default value of 5 Angstroms, otherwise Spinach would ignore long-range dipolar interactions.
Zero-field splitting	inter.coupling.matrix (as an interaction of the electron with itself) 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	inter.coupling.scalar
Inter-electron dipolar coupling	inter.coordinates if point dipolar approximation applies and electron coordinates are known (they will be converted into a dipolar interaction matrix internally), inter.coupling.matrix for dipolar coupling supplied as a matrix (make sure it's traceless), inter.coupling.eigs & inter.coupling.euler for dipolar interactions supplied as eigenvalues and Euler angles (make sure eigenvalues sum up to zero).

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 follows:

    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 sys.tols.dd_ncells field. The default value is to account for two images in each direction.

Liquid crystal order matrix

For partially oriented systems, the 3x3 Saupe order matrix may be supplied. Example:

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

At the 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 import

Interaction and spin system information may also be loaded from a Gaussian03/09 magnetic properties calculation log. The following options should be present in the route section of the Gaussian input file:

    #p nmr=(giao,spinspin) output=pickett

to enable the calculation of magnetic interaction parameters. Note also that for non-magnetic nuclei Gaussian assumes the nearest magnetic isotope. Import of the following properties is supported (if found in the log): coordinates, g-tensor, hyperfine coupling tensors, chemical shielding tensors, scalar couplings, quadrupolar couplings, magnetic susceptibility, spin-rotation tensors, SCF energy, Gibbs free energy.

Gaussian input processing is handled by two functions. The first one, gparse.m, reads every useful piece of information it can find in the log and returns that information as a Matlab data structure:

    props=gparse('..\molecules\strychnine.log');

Below is the list of parameters that gparse.m returns (all tensors are returned in the coordinate system corresponding to the “standard orientation” used by Gaussian03/09:

Gaussian parser output
props.inp_geom	Input orientation, XYZ, in Angstrom.
props.std_geom	Standard orientation, XYZ, in Angstrom.
props.natoms	Number of atoms in the molecule.
props.energy	SCF energy, in Hartrees.
props.method	Quantum chemical method used in evaluating the energy.
props.hfc.iso	Isotropic hyperfine couplings, in Gauss.
props.hfc.full.eigvals	Eigenvalues of the full hyperfine coupling tensor, in Gauss.
props.hfc.full.eigvecs	Eigenvectors of the full hyperfine coupling tensor.
props.hfc.full.matrix	Full hyperfine coupling tensors, in Gauss.
props.g_tensor.eigvals	Eigenvalues of the g-tensor, in Bohr magneton units.
props.g_tensor.eigvecs	Eigenvectors of the g-tensor.
props.g_tensor.matrix	Full g-tensor matrix, in Bohr magneton units.
props.cst	Chemical shielding tensors for each spin, in ppm.
props.nqi	Nuclear quadrupolar interaction tensors for each spin, Hz.
props.srt	Nuclear spin-rotation tensors for each spin, Hz.
props.chi	Magnetic susceptibility tensor, cubic Angstrom.
props.gibbs	Gibbs free energy.
props.j_coupling	Matrix of scalar couplings between each pair of spins, Hz.
props.symbols	Atomic symbols.
props.atomic_numbers	Atomic numbers.
props.filename	Name of the log file.

If a particular parameter is encountered multiple times in the Gaussian log, the last instance is returned. The function performs automatic symmetrisation of shielding tensors, hyperfine tensors and g-tensors. If unsymmetrized tensors are required, a second argument should be passed to the function:

    props=gparse('..\molecules\strychnine.log',options);

where options is a cell array of strings with the following strings permitted as elements: 'g_nosymm', 'cst_nosymm', 'hfc_nosymm'. The presence of any of these strings would cause the corresponding symmetrisation operation to be switched off.

After the log has been parsed, the second function, g2spinach.m, should be used to pick out the interactions that the user would like to import into the Spinach calculation:

    [sys,inter]=g2spinach(props,nuclei,references,options);

where props is the output of gparse.m and nuclei is a cell array of the following form:

    {{'H','1H'},{'N','15N'}...}

giving the list of elements and isotopes that should be included. All other nuclei will be ignored. The third argument is a vector of reference chemical shieldings in ppm. Gaussian03/09 returns absolute shielding (relative to the bare nucleus in vacuum) and it is therefore necessary to run a separate calculation for the reference substance (TMS for 1H and 13C) using the same method. The resulting reference chemical shielding values go into the references parameter.

Computed absolute isotropic shielding values for TMS in vacuum (ppm) for a few popular theory levels are:

     GIAO                  13C       1H 
     B3LYP/6-31G*        189.6621  32.1833 
     B3LYP/6-311+G(2d,p) 182.4485  31.8201 
     HF/6-31G*           199.9711  32.5957 
     HF/6-311+G(2d,p)    192.5828  32.0710

     CSGT                  13C       1H 
     B3LYP/6-31G*        188.5603  29.1952 
     B3LYP/6-311+G(2d,p) 182.1386  31.7788 
     HF/6-31G*           196.8670  29.5517 
     HF/6-311+G(2d,p)    192.5701  31.5989 

If you have other isotopes in the system, or a different reference substance, you will have to run Gaussian03 for that substance and use the resulting absolute shielding values. If the isotope list described above involves an electron, for example:

    {{'E','E'},{'H','1H'}...}

then ESR parameters are imported – chemical shielding and scalar couplings are ignored, but g-tensor and hyperfine couplings are included. A spin-1/2 electron shell is assumed in this case. A spin-1 (for example) electron may be specified by {'E','3E'}. When in doubt, it is usually best to look through props data structure manually and pick out the parameters that you wish to have in the simulation. The following options may be passed to g2spinach.m function:

Option	Units / Values	Comments
options.min_j	Hz	Scalar coupling threshold. J-couplings smaller than this value will be ignored.
options.min_hfc	Hz	Hyperfine coupling threshold. Hyperfine coupling tensors with 2-norm smaller than this value will be ignored.
options.purge	logical	If this switch is set, the nuclei with negligible hyperfine couplings are purged from the spin system.
options.no_xyz	logical	If this switch is set, the atomic coordinate information is ignored.

The following fields, if found in the original Gaussian log, are written by g2spinach.m function into the output parameters:

sys.isotopes	[nspins x 1] cell array of strings
inter.coordinates	[nspins x 1] cell array of 3-vectors, Angstrom. Not returned if there is an electron in the isotope list (because hyperfine interaction tensors already contain the through-space dipolar contribution to the electron-nuclear coupling).
inter.zeeman.matrix	[nspins x 1] cell array of 3x3 matrices, ppm for nuclei, Bohr magneton units for electrons.
inter.coupling.matrix	[nspins x nspins] cell array of 3x3 matrices, all in Hz.
inter.coupling.scalar	[nspins x nspins] cell array of scalars, all in Hz.
inter.spinrot.matrix	[nspins x 1] cell array of 3x3 matrices, all in Hz. Note that spin-rotation tensors are imported, but not supported by Spinach at the moment.

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 re-turns the corresponding Spinach data structures.

COSMOS import

Spin system information may also be read from a COSMOS *.coo file using the following syntax:

    [sys,inter]=cosmos2spinach(filename);

Cartesian coordinates and chemical shielding tensors are currently read from COSMOS log files. Chemical shielding tensors imported from COSMOS are not translated into chemical shift, this should be done manually.

Protein data import

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

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

would automatically form the necessary data structures. In this expression, pdb_file is a string specifying the PDB file (the first geometry is read if multiple geometries are present), bmrb_file is a string specifying the BMRB file and selection may be set to a list of numbers (in which case only spins with those numbers are loaded), to 'backbone' (in which case the protein backbone up to CB and HB is loaded) or to 'backbone-minimal' (in which case only the backbone, including CA and HA, is loaded). The function performs sophisticated post-processing, CSA assignment and J-coupling estimation. Particulars are not yet published – see the function code for more information about what exactly is going on.

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