Heteronuclear NMR simulations
Heteronuclear liquid-state NMR simulations in Spinach are set up in the same way as homonuclear simulations, but the pulse sequence parameters must specify multiple spin channels. This tutorial follows two examples from the Spinach example set: hsqc_strychnine.m, which simulates a phase-sensitive natural-abundance 13C HSQC spectrum of strychnine, and clip_hsqc_camphor.m, which simulates a natural abundance CLIP-HSQC spectrum of camphor.
HSQC spectrum of strychnine
Spinach examples are functions rather than scripts so that each run starts from a clean Matlab workspace. We must therefore open a new file and begin with with a function declaration:
function hsqc_strychnine()
Unlike the previous tutorial where we had specified the spin system manually, here we use one of the standard spin systems supplied with Spinach - the strychnine spin system:
[sys,inter]=strychnine({'1H','13C'});
The strychnine helper returns the sys and inter input structures for the selected isotopes. In this example, both proton and carbon-13 spins are imported. The strychnine data file contains proton shifts and couplings, carbon and nitrogen shifts and coordinates, and one-bond 13C-1H J-couplings. It does not provide 13C-13C J-couplings, and its source header states that it is suitable for natural-abundance 13C simulations.
The magnetic field is then set in Tesla:
sys.magnet=5.9;
The next lines request Spinach's greedy parallelisation option and set the proximity cut-off:
sys.enable={'greedy'}; sys.tols.prox_cutoff=4.0;
The greedy option is a parallel resource-allocation option: when Spinach starts a new parallel pool, every worker is allowed to use every CPU core. The proximity cut-off is in Angstrom and supplies a spatial-neighbour criterion to Spinach routines that require one.
HSQC basis set
The basis set is the same restricted spherical-tensor Liouville-space basis used by many liquid-state examples:
bas.formalism='sphten-liouv'; bas.approximation='IK-2'; bas.connectivity='scalar_couplings'; bas.space_level=1;
The sphten-liouv formalism is required by hsqc.m. The IK-2 approximation builds a restricted state space rather than the complete Liouville space, and scalar_couplings instructs the basis generator to use the scalar-coupling graph for connectivity. The space_level setting controls how far spatially local correlations are allowed into the basis.
HSQC sequence parameters
The sequence parameters in the example are:
parameters.J=140; parameters.sweep=[10000 3000]; parameters.offset=[4000 1000]; parameters.npoints=[128 128]; parameters.zerofill=[512 512]; parameters.spins={'13C','1H'}; parameters.decouple_f1={'1H'}; parameters.decouple_f2={'13C'}; parameters.axis_units='ppm';
The hsqc.m header specifies parameters.sweep, parameters.npoints, parameters.spins, parameters.decouple_f1, parameters.decouple_f2, and parameters.J. Here parameters.J is the working heteronuclear scalar coupling in Hz, used to set the transfer delay. The two entries in parameters.sweep are the F1 and F2 sweep widths in Hz, and the two entries in parameters.offset are the corresponding transmitter or receiver offsets in Hz. The point counts in parameters.npoints are the acquired F1 and F2 dimensions, and parameters.zerofill gives the Fourier transform sizes used during processing. The spin list {'13C','1H'} means that the indirect F1 dimension is carbon-13 and the directly detected F2 dimension is proton. The decouple_f1 field lists spins that receive midpoint 180-degree refocusing pulses during F1 evolution, and decouple_f2 lists spins decoupled during F2 acquisition. The plotted axes are requested in ppm.
HSQC spin-system construction
The Spinach spin-system object is created first:
spin_system=create(sys,inter);
The create function performs Spinach input processing and prints a report that should be checked for warnings. Natural-abundance carbon is handled by isotope dilution:
subsystems=dilute(spin_system,'13C');
The dilute function returns a cell array of spin systems. With the default tuple size, each returned subsystem contains a single instance of the specified dilute isotope and all other carbon-13 spins have been removed. This makes a natural-abundance 13C HSQC simulation practical: the final spectrum is obtained as a sum over all single-carbon isotopomers.
The example then preallocates the complex answer:
spectrum=zeros(parameters.zerofill(2),... parameters.zerofill(1),'like',1i);
The first dimension of the matrix corresponds to the direct F2 dimension, and the second dimension corresponds to the indirect F1 dimension. The use of 'like',1i requests a complex array.
HSQC isotopomer loop
Each isotopomer is processed independently:
parfor n=1:numel(subsystems)
The basis set is built inside the loop:
subsystem=basis(subsystems{n},bas);
This is necessary because dilute removes basis-set information when it changes the spin system. The HSQC simulation is then run in the liquid-state context with NMR assumptions:
fid=liquid(subsystem,@hsqc,parameters,'nmr');
The liquid context builds the isotropic Liouvillian, applies offsets and any other common sequence settings, and passes the resulting Hamiltonian, relaxation superoperator, and kinetics superoperator to hsqc.m. The HSQC sequence returns two components of the States quadrature signal:
fid.pos=apodisation(spin_system,fid.pos,{{'sqcos'},{'sqcos'}}); fid.neg=apodisation(spin_system,fid.neg,{{'sqcos'},{'sqcos'}});
Both components are multiplied by square-cosine windows in both dimensions. The direct dimension is transformed first:
f1_pos=fftshift(fft(fid.pos,parameters.zerofill(2),1),1); f1_neg=fftshift(fft(fid.neg,parameters.zerofill(2),1),1);
The States signal is then formed:
fid=f1_pos+conj(f1_neg);
and the indirect dimension is transformed and added to the accumulated spectrum:
spectrum=spectrum+fftshift(fft(fid,parameters.zerofill(1),2),2);
The loop is then closed:
end
Because each subsystem is independent, the loop is a natural use case for Matlab parfor.
HSQC plotting
The final stage opens a figure, scales it, and plots positive contours:
kfigure(); scale_figure([1.5 2.0]); plot_2d(spin_system,real(spectrum),parameters,... 20,[0.05 1.0 0.05 1.0],2,256,6,'positive');
The plot_2d arguments are the spin system, the real part of the spectrum, the experiment parameters, the number of contours, positive and negative contour elevation ranges, the non-linear contour-spacing curvature, the colour-map size, the colour-map curvature, and the contour sign selection. In this example the last argument requests positive contours only.
CLIP-HSQC spectrum of camphor
The clip_hsqc_camphor.m example also begins with a function declaration:
function clip_hsqc_camphor()
The spin system is imported from a Gaussian log file:
options.min_j=3.0; options.no_xyz=0;
[sys,inter]=g2spinach(gparse('../standard_systems/camphor.log'),...
{{'H','1H'},{'C','13C'}},[31.8 182.1],options);
The gparse function reads the Gaussian output file. The g2spinach function converts the parsed electronic-structure data into Spinach input structures. The particle list imports hydrogen atoms as 1H and carbon atoms as 13C. The reference vector [31.8 182.1] supplies the absolute shielding references used by the conversion to chemical shifts. The options.min_j=3.0 setting discards scalar couplings smaller than 3 Hz, and options.no_xyz=0 keeps the coordinate information.
The magnetic field is then specified:
sys.magnet=14.1;
The next block replaces the isotropic parts of the shielding tensors with experimental chemical shifts:
inter.zeeman.matrix=shift_iso(inter.zeeman.matrix,1:26,... [ 29.70 26.80 44.20 59.70 42.80 ... 218.1 49.70 17.80 17.30 7.80 ... 1.35 1.67 1.97 1.31 1.96 ... 0.85 0.85 0.85 0.99 0.99 ... 0.99 0.90 0.90 0.90 2.33 1.76]);
The shift_iso function preserves the anisotropic rank-1 and rank-2 parts of each tensor and replaces only the isotropic component with the supplied value. In this example, the source comment states that coordinates, shielding anisotropies, and J-couplings come from DFT, while isotropic chemical shifts come from experimental data.
CLIP-HSQC basis set and options
The basis set specification is:
bas.formalism='sphten-liouv'; bas.approximation='IK-2'; bas.connectivity='scalar_couplings'; bas.space_level=1;
The algorithmic options are:
sys.enable={'greedy'}; sys.tols.prox_cutoff=4.0;
As in the strychnine example, this requests Spinach's greedy parallelisation option and sets the proximity cut-off to 4 Angstrom.
CLIP-HSQC sequence parameters
The CLIP-HSQC example uses the following sequence parameters:
parameters.J=140; parameters.sweep=[8000 1500]; parameters.offset=[4000 1000]; parameters.npoints=[128 128]; parameters.zerofill=[512 512]; parameters.spins={'13C','1H'}; parameters.axis_units='ppm';
The clip_hsqc.m header requires parameters.sweep, parameters.npoints, parameters.spins, and parameters.J. As above, parameters.J is the working scalar coupling in Hz, parameters.sweep gives the F1 and F2 sweep widths, parameters.offset gives the corresponding offsets, parameters.npoints gives the acquired point counts, parameters.zerofill gives the Fourier transform sizes, parameters.spins={'13C','1H'} specifies carbon-13 in F1 and proton in F2, and parameters.axis_units='ppm' requests ppm axes.
CLIP-HSQC simulation
The spin-system construction and isotopomer generation steps are the same as in the HSQC example:
spin_system=create(sys,inter); subsystems=dilute(spin_system,'13C');
The answer is again preallocated as a complex matrix:
spectrum=zeros(parameters.zerofill(2),... parameters.zerofill(1),'like',1i);
The isotopomer loop starts with basis construction:
parfor n=1:numel(subsystems)
subsystem=basis(subsystems{n},bas);
The simulation call differs only in the pulse sequence function handle:
fid=liquid(subsystem,@clip_hsqc,parameters,'nmr');
The clip_hsqc.m sequence returns fid.pos and fid.neg as the two States quadrature components. The data processing matches the HSQC example:
fid.pos=apodisation(spin_system,fid.pos,{{'sqcos'},{'sqcos'}}); fid.neg=apodisation(spin_system,fid.neg,{{'sqcos'},{'sqcos'}}); f1_pos=fftshift(fft(fid.pos,parameters.zerofill(2),1),1); f1_neg=fftshift(fft(fid.neg,parameters.zerofill(2),1),1); fid=f1_pos+conj(f1_neg); spectrum=spectrum+fftshift(fft(fid,parameters.zerofill(1),2),2); end
The summation over the isotopomer loop produces the natural-abundance 13C CLIP-HSQC spectrum.
CLIP-HSQC plotting
The example plots negative contours:
kfigure(); scale_figure([1.5 2.0]); plot_2d(spin_system,real(spectrum),parameters,... 20,[0.05 0.5 0.05 0.5],2,256,6,'negative');
The contour-level vector differs from the strychnine example, and the last argument requests negative contours only.
Version 2.12, authors: Ilya Kuprov, Luke Edwards
