# Hyperfine shift module

Spinach includes a variety of functions that compute hyperfine shifts (including, contact, pseudocontact and dipolar components) and assist in their analysis. The formalism Spinach uses to compute them is described in our recent paper on the subject (http://dx.doi.org/10.1039/C4CP03106G). This section contains an overview of the related functions, located in kernel/utilities and experiments folders.

## Infrastructure functions

cgsppm2ang.m – converts magnetic susceptibility from the cgs-ppm (aka cm3/mol) units quoted by quantum chemistry packages into Å3 units required by Spinach pseudocontact shift functionality. Syntax:

    ang=cgsppm2ang(cgsppm)


Arrays of any dimension are supported.

fdhess.m – returns the finite-difference Hessian of a 3D array using a finite difference scheme with a user-specified number of stencil points. The dimensions of the 3D array are assumed to be ordered as [X Y Z]. Syntax:

    H=fdhess(A,npoints)


where the second parameter specifies the number of points in the finite difference stencil that should be used. The result is a 3x3 cell array of 3D matrices ordered in the following way:

    {d2A_dxdx  d2A_dxdy  d2A_dxdz
d2A_dydx  d2A_dydy  d2A_dydz
d2A_dzdx  d2A_dzdy  d2A_dzdz}


fdkup.m – returns a finite difference representation of the Kuprov operator:

    K[rho]=-(1/3)*Trace(Hessian[rho]*chi)


with the number of stencil points in the finite difference approximation specified by user. Syntax:

    K=fdkup(dims,chi,npoints)


where dims is a three-element vector specifying the dimensions of the 3D cube of data that the operator will be acting on, chi is the electron magnetic susceptibility tensor and npoints is the number of points in the finite difference stencil. The resulting operator is a sparse matrix designed to act on the vectorization of rho. The dimensions of rho are assumed to be ordered as [X Y Z].

fdlap.m – returns a finite-difference representation of the Laplacian for a 3D array with a user-specified finite difference stencil size. The resulting operator is a sparse matrix designed to act on the vectorization of the 3D array. The dimensions of the 3D array are assumed to be ordered as [X Y Z]. Syntax:

    L=fdlap(dims,npoints)


where dims is a three-element vector specifying the dimensions of 3D cube of data that the operator will be acting on, and npoints is the number of points in the finite difference stencil. Note: Dirichlet boundary conditions -- the resulting Laplacian should only be used for the inverse Laplacian operation.

fdmat.m – returns arbitrary-order central finite-difference differentiation matrix with unit spacing and sided finite-difference schemes of the same order at the edges. Syntax:

    D=fdmat(dim,npoints,order)


Parameters: dim - dimension of the column vector to be differentiated; npoints - number of points in the finite difference stencil; order - order of the derivative required.

interpmat.m – returns a matrix that acts on a stretched pseudocontact shift density cube and projects out the values of the PCS at the Cartesian coordinates given. Tricubic interpolation is used. Syntax:

    P=interpmat(cube_dims,ranges,xyz)


Parameters: cube_dims - PCS dimensions, ordered as [X Y Z]; ranges - Cartesian axis extents for the pseudocontact shift cube as [xmin xmax ymin ymax zmin zmax] in Angstroms; xyz - nuclear coordinates as [x y z] with multiple rows) at which PCS is to be evaluated, in Angstroms. Output: P - matrix projecting out PCS values at the given positions from the stretched PCS cube. Note: this function is a part of the PCS inverse problem solver module and should not normally be called directly by the user.

fdvec.m – performs arbitrary-order finite-difference differentiation of a user-supplied row or column vector. Uses central finite-difference stencils in the middle and sided stencils of the same order of accuracy on the sides. Syntax:

    dx=fdvec(x,npoints,order)


Parameters: x - column or row vector to be differentiated; npoints - number of points in the finite difference stencil; order - order of the derivative required.

fdweights.m – calculates finite difference weights for numerical derivatives (including order 0, which amounts to interpolation). Syntax:

    w=fdweights(target_point,grid_points,max_order)


Parameters: target_point - the point at which the derivative is required; grid_points - the points at which the function is given; max_order - maximum derivative order. The function returns finite difference coefficient array with the coefficients for the successive derivatives in rows.

points2mult.m - computes multipole moments from a set of point charges

   Ilm=points2mult(xyz,mxyz,rho,L,method)


Parameters: xyz - coordinates as [x y z] at which density rho is evaluated, in Angstroms, mxyz - coordinates of origin as [x y z], in Angstroms. L - array of ranks of spherical harmonics of the density L=[0 1 2], rho - column of the density at the points xyz, method - 'points' or 'grid'. 'Grid' integrates with corresponding d3r and 'points' just sum the density at points. Output: Ilm - predicted multipole moments of the density.

qform2sph.m - Returns the spherical expansion coefficients of RT*A*R/R^2

   [r0,r1,r2]=qform2sph(A)


Parameters: A - 3x3 matrix,Output: [r0,r1,r2] - coefficients for zero, first, and second rank spherical harmonics in the order -m:m

## Computing paramagnetic shifts

hfc2pcs.m – converts hyperfine coupling tensors and susceptibility tensors into pseudocontact shifts using Equation 10 from the paper by Charnock and Kuprov (http://dx.doi.org/10.1039/C4CP03106G). Syntax:

    sigma=hfc2pcs(A,chi,isotope)


Parameters: A - hyperfine coupling tensor, in MHz; chi - magnetic susceptibility tensor, in Angstrom3. Outputs: sigma - pseudocontact shift, in ppm.

hfc2pms.m – converts hyperfine coupling tensors and susceptibility tensors into paramagnetic shift tensors using Equation 10 from the paper by Charnock and Kuprov (http://dx.doi.org/10.1039/C4CP03106G). Syntax:

    sigma=hfc2pms(A,chi,isotope)


Parameters: A - hyperfine coupling tensor, in MHz; chi - magnetic susceptibility tensor, in Angstrom3. Outputs: sigma - dipolar shift tensor, in ppm. Note: both the susceptibility tensor and the hyperfine tensor must be submitted to this function with their isotropic parts preserved – do not make them traceless.

kpcs.m – computes the three-dimensional distribution of pseudocontact shift field by solving Kuprov equation for PCS. Syntax:

    [pcs_cube,pcs_vals]=kpcs(probden,chi,ranges,nxyz)


Parameters: probden - electron probability density cube; chi - electron magnetic susceptibility tensor in cubic Angstroms; ranges - Cartesian axis extents for the electron spin density cube as [xmin xmax ymin ymax zmin zmax] in angstroms; nxyz - nuclear coordinates as [x y z] with multiple rows) at which PCS is to be evaluated, in Angstroms. Output: pcs_vals - pseudocontact shift in ppm at each nucleus; pcs_cube - pseudocontact shift field on the same grid as the spin density supplied. Note: minimal three-point schemes are used for the finite difference operators. Increase stencil size if you have enough memory. Note: for further information on the equations and algorithms used in this function see http://dx.doi.org/10.1039/C4CP03106G.

ppcs.m – computes pseudocontact shift from a point electron centre at the nuclear coordinates supplied. Syntax:

    pcs_vals=ppcs(nxyz,mxyz,chi)


Parameters: chi - electron magnetic susceptibility tensor in cubic Angstroms; nxyz - nuclear coordinates as [x y z] with multiple rows, at which PCS is to be evaluated, in Angstroms; mxyz - paramagnetic centre coordinates as [x y z], in Angstroms; Output: pred_pcs - predicted pseudocontact shift (in ppm) at each of the nuclei.

lpcs.m - computes PCS from paramagnetic density (rho) outside the bounding sphere based on the multipole moments expansion

Parameters: nxyz - nuclear coordinates as [x y z] with multiple rows, at which PCS is to be evaluated, in Angstroms. mxyz - paramagnetic center coordinates as [x y z], in Angstroms. L - array of ranks of spherical harmonics of the density, Ilm - {[],[]} array of numbers corresponding to the integral Int(rho(r,theta,phi)Y*lm(theta,phi)r^ld^3r): for L=0 Ilm=N/2/sqrt(pi), for L=1 Ilm=[imag(I11) I10 real(I11)], for L=2 Ilm=[imag(I22) imag(I21) I20 real(I21) real(I22)], chi - 3x3 matrix or 5 values for susceptibility tensor in SI Angstrom^3. Output: theor_pcs - predicted pseudocontact shift (in ppm) at each of the nuclei.

## Visualization

volplot.m – volumetric plot function for scalar fields. Syntax:

    volplot(data_cube,axis_ranges)


where data_cube is a 3D cube of real data and axis ranges is a six-element vector giving axis extents as [xmin xmax ymin ymax zmin zmax].

## Fitting paramagnetic shifts

ipcs.m – solves the inverse problem for pseudocontact shift by recovering the source term in the Kuprov equation using Tikhonov + maximum entropy regularization procedure. Syntax:

    [source_cube,ranges,pred_pcs,diag_data]=ipcs(nxyz,expt_pcs,chi,npoints,lambda,margins)


Parameters: nxyz - nuclear coordinates as [x y z] with multiple rows at which PCS has been measured, in Angstroms; expt_pcs - pseudocontact shift in ppm at each nucleus; chi - electron magnetic susceptibility tensor, in units of Angstrom^3; npoints - number of points in each dimension of the source cube, a positive integer greater than 10; lambda - regularization parameters, the first element is the coefficient in front of the maximum entropy term and the second element is the coefficient in front of the Tikhonov term; margins - a six-element vector specifying margins to take around the bounding box of the nuclear coordinates supplied, to account for the possibility that the electron may be located on the periphery. Outputs: source_cube - source term cube with dimensions ordered as [X Y Z]; ranges - Cartesian axis extents for the source cube as [xmin xmax ymin ymax zmin zmax] in Angstroms; pred_pcs - pseudocontact shifts produced by the source cube returned in the first parameter; diag_data - the first element is the least squares error in ppm^2, the second element is the entropy penalty in the error functional, the third element is the tikhonov penalty in the error functional. Note: for further information on the equations and algorithms used in this function see http://dx.doi.org/10.1039/C4CP03106G

ippcs.m – fits the point electron model PCS to the experimental pseudocontact shift coordinates and values. Syntax:

    [exyz,chi,pred_pcs]=ippcs(nxyz,mguess,expt_pcs)


Parameters: nxyz - nuclear coordinates as [x y z] with multiple rows, at which PCS is measured, in Angstroms; mguess - initial guess for the unpaired electron coordinates as [x y z], in Angstroms; expt_pcs - pseudocontact shift in ppm at each nucleus. Outputs: exyz - optimized unpaired electron coordinates as [x y z], in Angstroms; chi - optimized magnetic susceptibility tensor in cubic, Angstroms; pred_pcs - predicted pseudocontact shift at each nucleus with the optimized mxyz and chi, ppm; Note: a good initial guess for the electron location is essential for a successful fit.

ilpcs.m - fits experimental PCS using the non-point approximation with multipole expansion of the density.

    mxyz,Ilm,pred_pcs]=ilpcs(nxyz,expt_pcs,L,mguess)


Parameters: nxyz - nuclear coordinates as [x y z] with multiple rows, at which PCS is to be evaluated, in Angstroms, expt_pcs - experimental pseudocontact shift (in ppm), L - row of ranks of spherical harmonics of the density. Output: mxyz - paramagnetic centre coordinates as [x y z], in Angstroms, Ilm - {[],[]} cell array of numbers corresponding to the multipole moments: Int(rho(r,theta,phi)Ylm(theta,phi)r^(l+2)sin(theta)d^3r): for L=0 Ilm=N/2/sqrt(pi), for L=1 Ilm=[real(I11) I10 imag(I11)], for L=2 Ilm=[real(I22) real(I21) I20 imag(I21) imag(I22)], chi - susceptibility tensor in SI Angstrom^3, pred_pcs - predicted pseudocontact shift (in ppm) at each of the nuclei.

Version 2.1, authors: Ilya Kuprov