Difference between revisions of "Spin system specification"

From Spinach Documentation Wiki
Jump to: navigation, search
(→‎Spin-spin interactions)
(Document the inter.modes bosonic mode input block (Spinach 2.12, PR 195))
 
(40 intermediate revisions by 2 users not shown)
Line 3: Line 3:
 
     [[Appendix D: kernel data structures|spin_system]]=[[create.m|create]]([[sys]],[[inter]]);
 
     [[Appendix D: kernel data structures|spin_system]]=[[create.m|create]]([[sys]],[[inter]]);
  
−
Follow the instructions below to prepare [[sys]] and [[inter]] arguments. A somewhat easier alternative is to use the [http://spindynamics.org/Spinach-GUI.php 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.
+
Follow the instructions below to prepare [[sys]] and [[inter]] arguments.
  
 
==Isotopes and text labels==
 
==Isotopes and text labels==
Line 12: Line 12:
  
 
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.
 
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'.
+
''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'.
 +
 
 +
Cavity modes, phonon modes, and transmons are supported as bosonic particles: use 'C#', 'V#', and 'T#' notation, for example 'C10' requests a cavity mode with 10 population number levels. Their frequencies, anharmonicities, damping, and couplings are specified in the [[#Bosonic modes|bosonic modes]] block described below.
  
 
Optionally, specify a label for each spin by giving a list of strings, for example:
 
Optionally, specify a label for each spin by giving a list of strings, for example:
Line 19: Line 21:
 
     [[sys]].labels={'CA','CB','HB2','HB3'};
 
     [[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.
+
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.
 +
 
 +
To make isotope replacements in existing [[sys]] and [[inter]] data structures, use [[isoswap.m]] function.
  
 
==Zeeman interactions==
 
==Zeeman interactions==
Line 38: Line 42:
 
|[[inter]].zeeman.eigs, [[inter]].zeeman.euler
 
|[[inter]].zeeman.eigs, [[inter]].zeeman.euler
 
|[1 × nspins] cell arrays of [1 × 3] row vectors
 
|[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.
+
|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
 
|[[inter]].zeeman.matrix
 
|[1 × nspins] cell array of [3 × 3] matrices
 
|[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.
+
|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
 
|[[inter]].zeeman.scalar
 
|[1 × nspins] cell array of real numbers
 
|[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.
+
|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.
 
|}
 
|}
  
Line 63: Line 67:
  
 
==Spin-spin interactions==
 
==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.
+
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.
  
 
{| class="wikitable"
 
{| class="wikitable"
Line 84: Line 88:
 
|[[inter]].coordinates
 
|[[inter]].coordinates
 
|[nspins × 1] cell array of [1 × 3] row vectors
 
|[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.
+
|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
 
|[[inter]].giant.coeff
−
|A cell array with a three-element vector followed by a five-element vector, followed by a seven-element vector and so on
+
|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 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.
+
|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
 
|[[inter]].giant.euler
−
|[1 × nspins] cell array of [1 × 3] row vectors
+
|[1 × nspins] cell array of cell arrays of [1 × 3] row vectors
−
|Euler angles (in radians) specifying the molecular frame orientation of the interaction specified via the irreducible spherical tensor coefficients given in [[inter]].giant.coeff field.  
+
|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.  
 
|}
 
|}
  
−
Example:
+
Examples:
  
 
     [[inter]].coupling.scalar={0 50; 0 0};
 
     [[inter]].coupling.scalar={0 50; 0 0};
Line 102: Line 106:
 
     [[inter]].coordinates={[0.0 0.0 0.0] ...
 
     [[inter]].coordinates={[0.0 0.0 0.0] ...
 
                         [0.0 0.0 1.5]};
 
                         [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.
 
Spin-spin interactions may be specified in a variety of equivalent ways. The table below provides suggestions on specifying all common interactions.
Line 134: Line 141:
 
!Inter-electron dipolar coupling
 
!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 [[Appendix_A:_internal_tolerances|tolerance]]. The default is 100 Angstrom.
 
|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 [[Appendix_A:_internal_tolerances|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.
 +
|}
 +
 +
==Bosonic modes==
 +
 +
''Spinach'' treats electromagnetic cavity modes, phonon modes, and transmons as particles alongside spins: they are declared in [[sys]].isotopes with a label that carries the Fock space truncation level, and their frequencies, anharmonicities, damping, and couplings are declared in the [[inter]].modes structure. This is the input side of cavity and circuit QED, optomechanics, and spin-phonon dynamics; the corresponding simulation frames are described on the [[assume.m]] page and the pulse sequence interface on the [[device.m]] page.
 +
 +
{| class="wikitable"
 +
!colspan="3"|Bosonic particle labels
 +
|-
 +
!Label!!Particle!!Typical use
 +
|-
 +
|'C#'
 +
|electromagnetic cavity mode with # population levels
 +
|cavity and circuit QED, optical modes
 +
|-
 +
|'V#'
 +
|phonon or vibrational mode with # population levels
 +
|spin-phonon coupling, optomechanics
 +
|-
 +
|'T#'
 +
|transmon with # energy levels
 +
|superconducting qubits
 +
|}
 +
 +
For example:
 +
 +
    [[sys]].isotopes={'E','C5','V7','T3'};
 +
 +
declares an electron, a five-level cavity mode, a seven-level phonon mode, and a three-level transmon. At least three levels are required for each mode. The truncation level is a convergence parameter and it is the user's responsibility to check it: a cavity expected to hold a coherent state of amplitude alpha needs several levels above |alpha|^2, and a transmon used as a qubit usually needs three levels so that leakage into the second excited state is visible.
 +
 +
Bosonic particles have zero magnetogyric ratio: they acquire no Zeeman interaction from [[sys]].magnet, and [[create.m]] refuses [[inter]].zeeman, [[inter]].coupling, and [[inter]].coordinates entries for them. Everything they do is declared in [[inter]].modes, whose fields are all cell arrays with one element (or one row and column) per particle, empty for particles to which the field does not apply - the same convention that the Zeeman and coupling specifications use. All frequencies and couplings are specified in Hz.
 +
 +
===Mode energies===
 +
 +
{| class="wikitable"
 +
!Variable name!!Variable type!!Content
 +
|-
 +
|[[inter]].modes.frqs
 +
|[1 × nspins] cell array of real numbers
 +
|Mode frequency in Hz. Mandatory whenever [[inter]].modes is present, and mandatory for every bosonic particle. Negative values are only legal for transmons, or when a carrier is declared for that mode, in which case the number is a detuning. Zero means the mode is already in its own resonant rotating frame.
 +
|-
 +
|[[inter]].modes.carriers
 +
|[1 × nspins] cell array of positive numbers
 +
|Laboratory frequency in Hz of the rotating frame in which the corresponding element of frqs is specified. Declaring it makes the frequency a detuning, permits negative detunings, supplies the physical frequency to the thermal occupation calculation, and suppresses the laboratory-frequency notice for modes that are not coupled to any spin.
 +
|-
 +
|[[inter]].modes.anharms
 +
|[1 × nspins] cell array of real numbers
 +
|Duffing anharmonicity in Hz. Physical transmons have negative anharmonicity; a positive value on a transmon produces a warning. [[ejec2duffing.m]] converts Josephson and charging energies into a frequency and an anharmonicity.
 +
|}
 +
 +
===Damping and dephasing===
 +
 +
{| class="wikitable"
 +
!Variable name!!Variable type!!Content
 +
|-
 +
|[[inter]].modes.lifetimes
 +
|[1 × nspins] cell array of positive numbers
 +
|Amplitude damping lifetime in seconds; the damping rate is kappa=1/T1.
 +
|-
 +
|[[inter]].modes.linewidths
 +
|[1 × nspins] cell array of positive numbers
 +
|Amplitude damping linewidth in Hz; the damping rate is kappa=2*pi*linewidth.
 +
|-
 +
|[[inter]].modes.qfactors
 +
|[1 × nspins] cell array of positive numbers
 +
|Quality factor; the damping rate is kappa=|omega|/Q. A non-zero mode frequency is required, because on a zero-frequency mode this would silently produce zero damping; use lifetimes or linewidths there instead.
 +
|-
 +
|[[inter]].modes.t2_times
 +
|[1 × nspins] cell array of positive numbers
 +
|Total coherence time in seconds. The pure dephasing rate is extracted as 1/T2-kappa*(1+2*nbar)/2 at the declared temperature; a T2 longer than the damping limit allows is an error that reports the bound.
 +
|}
 +
 +
The three damping specifications are alternatives: declaring more than one of them for the same mode is an error. [[inter]].temperature is mandatory when any mode is damped, because the thermal occupation of the mode bath must be known; the room temperature default is refused with an error message naming the occupation number that it would have implied. The dissipators themselves are built by [[rlx_modes.m]] and require a Liouville space formalism.
 +
 +
===Couplings===
 +
 +
Four channels describe direct interactions between modes and between modes and spins. Each is an [nspins × nspins] cell array of real numbers in Hz.
 +
 +
{| class="wikitable"
 +
!Variable name!!Operator!!Content and placement
 +
|-
 +
|[[inter]].modes.exchange
 +
|a<sub>k</sub><sup>+</sup>a<sub>l</sub> + a<sub>k</sub>a<sub>l</sub><sup>+</sup>, or L<sub>+</sub>a + L<sub>-</sub>a<sup>+</sup>
 +
|Beam-splitter coupling between two modes, or Jaynes-Cummings coupling between a spin and a mode. Upper triangle only. At least one of the two particles must be a mode: spin-spin couplings belong in [[inter]].coupling.
 +
|-
 +
|[[inter]].modes.kerr
 +
|n<sub>k</sub>n<sub>l</sub>
 +
|Cross-Kerr coupling between two modes. Upper triangle only. Both particles must be bosonic modes.
 +
|-
 +
|[[inter]].modes.longitudinal
 +
|L<sub>Z</sub>(a+a<sup>+</sup>)/sqrt(2), or n(a+a<sup>+</sup>)/sqrt(2)
 +
|Coupling of a mode quadrature to a spin projection (longitudinal spin-phonon coupling) or to the population of another mode (optomechanical radiation pressure). One entry per pair, in either triangle: the row index carries the spin projection or population operator and the column index carries the quadrature. Declaring both {n,k} and {k,n} is an ambiguity error, as is a diagonal entry.
 +
|-
 +
|[[inter]].modes.dispersive
 +
|L<sub>Z</sub>n
 +
|Dispersive coupling of a spin projection to a mode population, as used for circuit QED readout. Upper triangle only. Exactly one spin and one mode; number-number couplings between two modes belong in the Kerr channel.
 +
|}
 +
 +
Couplings whose magnitude falls below [[Appendix A: internal tolerances|inter_cutoff]] are dropped, and the drop is reported to the console by channel and particle pair.
 +
 +
===Modulation of the spin Hamiltonian by mode coordinates===
 +
 +
Two further fields describe the modulation of the ''spin'' Hamiltonian by mode displacement, which is the mechanism behind spin-phonon relaxation and strain-driven spin dynamics. Both are [nspins × nspins] cell arrays indexed by mode pairs on or above the diagonal: a diagonal entry {k,k} holds derivatives with respect to the coordinate of mode k, and an off-diagonal entry {k1,k2} holds the mixed second derivative of a two-phonon Raman term. Each entry is a cell array of at most two elements, indexed by derivative order.
 +
 +
{| class="wikitable"
 +
!Variable name!!Variable type!!Content
 +
|-
 +
|[[inter]].modes.coupling_mod
 +
|nested cell arrays of [3 × 3] matrices
 +
|Derivatives of spin-spin coupling tensors in Hz with respect to dimensionless mode coordinates. The block for each derivative order is an [nspins × nspins] cell array indexed by spin pairs. Same-spin blocks must be traceless and symmetric; they are expanded into second-rank irreducible spherical tensors.
 +
|-
 +
|[[inter]].modes.zeeman_mod
 +
|nested cell arrays of [1 × 3] vectors
 +
|Derivatives of effective local fields in Hz with respect to dimensionless mode coordinates. The block for each derivative order is a row cell array with one element per particle. Column cell arrays are rejected.
 
|}
 
|}
 +
 +
First derivatives describe a single mode and therefore require a diagonal entry. The dimensionless mode coordinate is (a+a<sup>+</sup>)/sqrt(2); [[cart2mode.m]] converts derivatives with respect to Cartesian atomic displacements, as produced by electronic structure packages, into this convention, and [[mev2hz.m]], [[icm2hz.m]], and [[kelvin2hz.m]] convert the energy units in which such data usually arrives.
 +
 +
===Examples===
 +
 +
A spin resonant with a cavity mode, coupled by a Jaynes-Cummings interaction:
 +
 +
    % An electron and a five-level cavity mode
 +
    [[sys]].magnet=0.33;
 +
    [[sys]].isotopes={'E','C5'};
 +
 +
    % Cavity resonant with the electron Larmor frequency
 +
    [[inter]].modes.frqs={[] -sys.magnet*spin('E')/(2*pi)};
 +
 +
    % Jaynes-Cummings coupling of 2.828 MHz
 +
    [[inter]].modes.exchange=cell(2,2);
 +
    [[inter]].modes.exchange{1,2}=2.828e6;
 +
 +
A leaky microwave resonator at 50 mK, with an amplitude damping lifetime and a total coherence time:
 +
 +
    % Microwave resonator with five Fock levels
 +
    [[sys]].isotopes={'C5'};
 +
    [[inter]].modes.frqs={6.02e9};
 +
    [[inter]].modes.lifetimes={10e-9};
 +
    [[inter]].modes.t2_times={9.9e-9};
 +
    [[inter]].temperature=0.050;
 +
 +
A phonon mode coupled longitudinally to an electron, the input of a spin-phonon dephasing calculation:
 +
 +
    % An electron and a seven-level phonon mode
 +
    [[sys]].isotopes={'E','V7'};
 +
    [[inter]].modes.frqs={[] 20e6};
 +
 +
    % Longitudinal coupling: row 1 carries the spin
 +
    % projection, column 2 carries the quadrature
 +
    [[inter]].modes.longitudinal=cell(2,2);
 +
    [[inter]].modes.longitudinal{1,2}=4e6*sqrt(2);
 +
 +
When the mode block is ingested, [[create.m]] prints three console tables: mode parameters, mode couplings by channel, and the norms of the spin Hamiltonian modulation blocks. Reading those tables is the fastest way to confirm that the input was understood as intended.
  
 
==Periodic boundary conditions==
 
==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:
+
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]};
 
     [[inter]].pbc={[0 0 10]};
Line 151: Line 314:
  
 
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 [[Appendix_A:_internal_tolerances|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 [[Appendix_A:_internal_tolerances|tolerance]] settings.
 
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 [[Appendix_A:_internal_tolerances|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 [[Appendix_A:_internal_tolerances|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==
 
==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:
+
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]);
+
     [[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.
 
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==
 
==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:
+
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
 
     % Parse a DFT calculation log
Line 169: Line 350:
 
     <nowiki>[</nowiki>[[sys]],[[inter]]]=[[g2spinach.m|g2spinach]](props,{{'C','13C'},{'N','15N'}},[182.1 264.5],[]);
 
     <nowiki>[</nowiki>[[sys]],[[inter]]]=[[g2spinach.m|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.
+
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.
  
−
==SIMPSON import==
+
'''Some examples of Gaussian inputs:''' [[:File:gaussian_a.txt|gaussian_a.txt]], [[:File:gaussian_b.txt|gaussian_b.txt]], [[:File:gaussian_c.txt|gaussian_c.txt]] - note that the extra printing flag (#p) should always be present.
−
Spin system information may also be read from the spinsys{} field of a SIMPSON *.in file using the following syntax:
 
  
−
    <nowiki>[</nowiki>[[sys]],[[inter]]]=[[s2spinach.m|s2spinach]](filename);
+
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.
  
−
The function absorbs isotope information, all Zeeman tensors and all coupling tensors, and returns ''Spinach'' [[sys]] and [[inter]] data structures with that information.
+
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.
  
 
==Protein and nucleic acid import==
 
==Protein and nucleic acid import==
Line 239: Line 419:
 
On suitable hardware (solid state disks, NVidia Tesla cards), these options can significantly improve performance. See the [[Appendix_F:_expert-level_options|expert level options]] page for further information.
 
On suitable hardware (solid state disks, NVidia Tesla cards), these options can significantly improve performance. See the [[Appendix_F:_expert-level_options|expert level options]] page for further information.
  
−
 
+
''Version 2.8, authors: [[Ilya Kuprov]], [[Luke Edwards]], [[Hannah Hogben]]''
−
''Version 1.9, authors: [[Ilya Kuprov]], [[Luke Edwards]], [[Hannah Hogben]]''
 

Latest revision as of 05:43, 6 August 2026

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

Cavity modes, phonon modes, and transmons are supported as bosonic particles: use 'C#', 'V#', and 'T#' notation, for example 'C10' requests a cavity mode with 10 population number levels. Their frequencies, anharmonicities, damping, and couplings are specified in the bosonic modes block described below.

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.

To make isotope replacements in existing sys and inter data structures, use isoswap.m function.

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.

Bosonic modes

Spinach treats electromagnetic cavity modes, phonon modes, and transmons as particles alongside spins: they are declared in sys.isotopes with a label that carries the Fock space truncation level, and their frequencies, anharmonicities, damping, and couplings are declared in the inter.modes structure. This is the input side of cavity and circuit QED, optomechanics, and spin-phonon dynamics; the corresponding simulation frames are described on the assume.m page and the pulse sequence interface on the device.m page.

Bosonic particle labels
Label Particle Typical use
'C#' electromagnetic cavity mode with # population levels cavity and circuit QED, optical modes
'V#' phonon or vibrational mode with # population levels spin-phonon coupling, optomechanics
'T#' transmon with # energy levels superconducting qubits

For example: 

    sys.isotopes={'E','C5','V7','T3'};

declares an electron, a five-level cavity mode, a seven-level phonon mode, and a three-level transmon. At least three levels are required for each mode. The truncation level is a convergence parameter and it is the user's responsibility to check it: a cavity expected to hold a coherent state of amplitude alpha needs several levels above |alpha|^2, and a transmon used as a qubit usually needs three levels so that leakage into the second excited state is visible.

Bosonic particles have zero magnetogyric ratio: they acquire no Zeeman interaction from sys.magnet, and create.m refuses inter.zeeman, inter.coupling, and inter.coordinates entries for them. Everything they do is declared in inter.modes, whose fields are all cell arrays with one element (or one row and column) per particle, empty for particles to which the field does not apply - the same convention that the Zeeman and coupling specifications use. All frequencies and couplings are specified in Hz.

Mode energies

Variable name Variable type Content
inter.modes.frqs [1 × nspins] cell array of real numbers Mode frequency in Hz. Mandatory whenever inter.modes is present, and mandatory for every bosonic particle. Negative values are only legal for transmons, or when a carrier is declared for that mode, in which case the number is a detuning. Zero means the mode is already in its own resonant rotating frame.
inter.modes.carriers [1 × nspins] cell array of positive numbers Laboratory frequency in Hz of the rotating frame in which the corresponding element of frqs is specified. Declaring it makes the frequency a detuning, permits negative detunings, supplies the physical frequency to the thermal occupation calculation, and suppresses the laboratory-frequency notice for modes that are not coupled to any spin.
inter.modes.anharms [1 × nspins] cell array of real numbers Duffing anharmonicity in Hz. Physical transmons have negative anharmonicity; a positive value on a transmon produces a warning. ejec2duffing.m converts Josephson and charging energies into a frequency and an anharmonicity.

Damping and dephasing

Variable name Variable type Content
inter.modes.lifetimes [1 × nspins] cell array of positive numbers Amplitude damping lifetime in seconds; the damping rate is kappa=1/T1.
inter.modes.linewidths [1 × nspins] cell array of positive numbers Amplitude damping linewidth in Hz; the damping rate is kappa=2*pi*linewidth.
inter.modes.qfactors [1 × nspins] cell array of positive numbers omega|/Q. A non-zero mode frequency is required, because on a zero-frequency mode this would silently produce zero damping; use lifetimes or linewidths there instead.
inter.modes.t2_times [1 × nspins] cell array of positive numbers Total coherence time in seconds. The pure dephasing rate is extracted as 1/T2-kappa*(1+2*nbar)/2 at the declared temperature; a T2 longer than the damping limit allows is an error that reports the bound.

The three damping specifications are alternatives: declaring more than one of them for the same mode is an error. inter.temperature is mandatory when any mode is damped, because the thermal occupation of the mode bath must be known; the room temperature default is refused with an error message naming the occupation number that it would have implied. The dissipators themselves are built by rlx_modes.m and require a Liouville space formalism.

Couplings

Four channels describe direct interactions between modes and between modes and spins. Each is an [nspins × nspins] cell array of real numbers in Hz.

Variable name Operator Content and placement
inter.modes.exchange ak+al + akal+, or L+a + L-a+ Beam-splitter coupling between two modes, or Jaynes-Cummings coupling between a spin and a mode. Upper triangle only. At least one of the two particles must be a mode: spin-spin couplings belong in inter.coupling.
inter.modes.kerr nknl Cross-Kerr coupling between two modes. Upper triangle only. Both particles must be bosonic modes.
inter.modes.longitudinal LZ(a+a+)/sqrt(2), or n(a+a+)/sqrt(2) Coupling of a mode quadrature to a spin projection (longitudinal spin-phonon coupling) or to the population of another mode (optomechanical radiation pressure). One entry per pair, in either triangle: the row index carries the spin projection or population operator and the column index carries the quadrature. Declaring both {n,k} and {k,n} is an ambiguity error, as is a diagonal entry.
inter.modes.dispersive LZn Dispersive coupling of a spin projection to a mode population, as used for circuit QED readout. Upper triangle only. Exactly one spin and one mode; number-number couplings between two modes belong in the Kerr channel.

Couplings whose magnitude falls below inter_cutoff are dropped, and the drop is reported to the console by channel and particle pair.

Modulation of the spin Hamiltonian by mode coordinates

Two further fields describe the modulation of the spin Hamiltonian by mode displacement, which is the mechanism behind spin-phonon relaxation and strain-driven spin dynamics. Both are [nspins × nspins] cell arrays indexed by mode pairs on or above the diagonal: a diagonal entry {k,k} holds derivatives with respect to the coordinate of mode k, and an off-diagonal entry {k1,k2} holds the mixed second derivative of a two-phonon Raman term. Each entry is a cell array of at most two elements, indexed by derivative order.

Variable name Variable type Content
inter.modes.coupling_mod nested cell arrays of [3 × 3] matrices Derivatives of spin-spin coupling tensors in Hz with respect to dimensionless mode coordinates. The block for each derivative order is an [nspins × nspins] cell array indexed by spin pairs. Same-spin blocks must be traceless and symmetric; they are expanded into second-rank irreducible spherical tensors.
inter.modes.zeeman_mod nested cell arrays of [1 × 3] vectors Derivatives of effective local fields in Hz with respect to dimensionless mode coordinates. The block for each derivative order is a row cell array with one element per particle. Column cell arrays are rejected.

First derivatives describe a single mode and therefore require a diagonal entry. The dimensionless mode coordinate is (a+a+)/sqrt(2); cart2mode.m converts derivatives with respect to Cartesian atomic displacements, as produced by electronic structure packages, into this convention, and mev2hz.m, icm2hz.m, and kelvin2hz.m convert the energy units in which such data usually arrives.

Examples

A spin resonant with a cavity mode, coupled by a Jaynes-Cummings interaction: 

    % An electron and a five-level cavity mode
    sys.magnet=0.33;
    sys.isotopes={'E','C5'};

    % Cavity resonant with the electron Larmor frequency
    inter.modes.frqs={[] -sys.magnet*spin('E')/(2*pi)};

    % Jaynes-Cummings coupling of 2.828 MHz
    inter.modes.exchange=cell(2,2);
    inter.modes.exchange{1,2}=2.828e6;

A leaky microwave resonator at 50 mK, with an amplitude damping lifetime and a total coherence time: 

    % Microwave resonator with five Fock levels
    sys.isotopes={'C5'};
    inter.modes.frqs={6.02e9};
    inter.modes.lifetimes={10e-9};
    inter.modes.t2_times={9.9e-9};
    inter.temperature=0.050;

A phonon mode coupled longitudinally to an electron, the input of a spin-phonon dephasing calculation: 

    % An electron and a seven-level phonon mode
    sys.isotopes={'E','V7'};
    inter.modes.frqs={[] 20e6};

    % Longitudinal coupling: row 1 carries the spin
    % projection, column 2 carries the quadrature
    inter.modes.longitudinal=cell(2,2);
    inter.modes.longitudinal{1,2}=4e6*sqrt(2);

When the mode block is ingested, create.m prints three console tables: mode parameters, mode couplings by channel, and the norms of the spin Hamiltonian modulation blocks. Reading those tables is the fastest way to confirm that the input was understood as intended.

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.

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

Retrieved from "https://spindynamics.org/wiki/index.php?title=Spin_system_specification&oldid=9888"