Simple liquid state NMR simulations

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

As a simple example that is both sufficiently easy to get you started and sufficiently difficult to require Spinach, consider the simulation of a magnitude-mode COSY45 spectrum of rotenone -- a system with 22 spins and an irregular coupling pattern. This example is available in the example set supplied with Spinach.

Chemical structure of rotenone.

The first thing Spinach requires is the function declaration:

    function cosy45_rotenone()


This is not strictly necessary, but a good practice because this guarantees that Matlab starts the simulation with a clean background where no previously assigned variables exist.

The second stage is to specify the isotopes, 22 protons in this case:

    sys.isotopes={'1H','1H','1H','1H','1H','1H','1H','1H','1H',...
'1H','1H','1H','1H','1H','1H','1H','1H','1H',...
'1H','1H','1H','1H'};


See the spin system specification section for details on how to specify spin systems. The next step is to specify the magnet field:

    sys.magnet=5.9;


then chemical shifts, in ppm for all protons:

    inter.zeeman.scalar={6.72 6.40 4.13 4.56 4.89 6.46 7.79 3.79 2.91...
3.27 5.19 4.89 5.03 1.72 1.72 1.72 3.72 3.72...
3.72 3.76 3.76 3.76};


then all scalar couplings, in Hz:

    inter.coupling.scalar{3,4}=12.1;
inter.coupling.scalar{4,5}=3.1;
inter.coupling.scalar{3,5}=1.0;
inter.coupling.scalar{3,8}=1.0;
inter.coupling.scalar{1,8}=1.0;
inter.coupling.scalar{6,7}=8.6;
inter.coupling.scalar{5,8}=4.1;
inter.coupling.scalar{7,9}=0.7;
inter.coupling.scalar{7,10}=0.7;
inter.coupling.scalar{9,10}=15.8;
inter.coupling.scalar{10,11}=9.8;
inter.coupling.scalar{9,11}=8.1;
inter.coupling.scalar{13,14}=1.5;
inter.coupling.scalar{12,14}=0.9;
inter.coupling.scalar{22,22}=0;


where the last line is necessary to tell Matlab that the array is 22 by 22 and all other elements are empty or zero. This completes the spin system specifications for situations where accurate relaxation theory treatment is not required.

The next stage is basis set specification. The complete basis set for a 22-spin system is too large and we must therefore rely on the restricted state space approximation. See the basis set specification section for further details on basis set selection. Here we will be using IK-2 basis set with Liouville space formalism and no spatial proximity analysis because coordinates are not supplied:

    bas.formalism='sphten-liouv';
bas.approximation='IK-2';
bas.space_level=1;
bas.connectivity='scalar_couplings';


The three methyl groups contain magnetically equivalent protons and this symmetry may optionally be used to reduce the calculation time:

    bas.sym_group={'S3','S3','S3'};
bas.sym_spins={[14 15 16],[17 18 19],[20 21 22]};


This completes the basis set specification. The next stage is to specify the pulse sequence parameters. The full list of what Spinach pulse sequences could require is located here. The specific parameters required by the COSY sequence are in the corresponding section of this manual as well as in the function header of the COSY sequence file. The parameters required by COSY are:

    parameters.angle=pi/4;
parameters.offset=1200;
parameters.sweep=2000;
parameters.npoints=[512 512];
parameters.zerofill=[2048 2048];
parameters.spins={'1H'};
parameters.axis_units='ppm';


where the names are self-explanatory. This completes the specification of the spin system, of the basis set and of the experiment parameters. The next stage is to give all that information to Spinach. This is accomplished by running the two housekeeping functions:

    spin_system=create(sys,inter);
spin_system=basis(spin_system,bas);


Both print copious output to the console which you should always read carefully because it might contain warning messages. The next stage is simulation, which we will carry out in liquid state (hence the liquid context function) with the kernel assumptions set to 'nmr', indicating common high-field NMR spectroscopy:

    fid=liquid(spin_system,@cosy,parameters,'nmr');


The simulation returns the two-dimensional free induction decay that should undergo apodization (cosine bell in both dimensions is a good choice here):

    fid=apodization(fid,'cosbell-2d');


and Fourier transform:

    spectrum=fftshift(fft2(fid,parameters.zerofill(2),parameters.zerofill(1)));


The streaks in the magnitude-mode spectrum may be reduced somewhat by the destreaking function:

    spectrum=destreak(spectrum);


Finally, the spectrum is plotted:

    plot_2d(spin_system,abs(spectrum),parameters,20,[0.0025 0.05 0.0025 0.05],2,256,6,'positive');


The whole simulation should take less than a minute on any modern PC. Note that Matlab auto-starts the parallelization engine when it runs for the first time, that stage only happens once per Matlab session.

COSY45 simulation for rotenone.

Version 1.9, authors: Ilya Kuprov, Matthew Krzystyniak