Error using alpha

Topics related to Spinach package
andrea_simion
Posts: 11
Joined: Sun Aug 15, 2021 9:42 pm

Error using alpha

Post by andrea_simion »

Hello,

After I wrote the spin system specifications, I created the spin system, I specified the basis, and it works. But when I try to obtain the hamiltonian or the relaxation superoperators (see code below), I obtain the following error:

Error using alpha
Too many output arguments.

Where is the problem?

Code: Select all


sys.magnet=11.74;
parameters.rate=20000;
parameters.axis=[1 1 1];

% chennels: channel 1: 1H, channel 2: 13C
parameters.spins={'1H','13C'};
  
% isotopes 
sys.isotopes={'13C','1H','1H','1H','1H','1H','1H','1H'};

% atomic_coords
inter.coordinates={[-3.153 1.989 -3.728];
                   [-3.4   3.047 -3.831];
                   [-4.104 0.083 -4.015];
                   [-5.146 1.294 -3.299];
                   [-4.76  1.363 -5.014];
                   [-1.743 0.731 -4.624];
                   [-2.442 1.836 -5.68];
                   [-1.3   2.363 -4.562]};

% Isotropic chemical shifts (in ppm)
inter.zeeman.scalar={0 0 -2.6 -2.6 -2.6 5.3 5.3 5.3};

% Full chemical shift tensor (in ppm)
inter.zeeman.matrix={[0 0 0; 0 0 0; 0 0 0],...
                     [-3.64 0 0; 0 -1.56 0; 0 0 5.2],...
                     [-0.44 0 0; 0 -9.16 0; 0 0 1.8],...
                     [-5.42 0 0; 0 -4.48 0; 0 0 2.1],...
                     [-4.64 0 0; 0 -4.23 0; 0 0 1.1],...
                     [-1.34 0 0; 0 -1.46 0; 0 0 18.7],...
                     [-5.02 0 0; 0 -1.58 0; 0 0 22.5],...
                     [-2.68 0 0; 0 -1.23 0; 0 0 19.8]};

% Euler angles (degrees expressed in radians)
inter.zeeman.euler={[0          0           0];
                    [0          20*pi/180   10*pi/180];
                    [90*pi/180  80*pi/180   60*pi/180];
                    [120*pi/180 130*pi/180  50*pi/180];
                    [50*pi/180  72*pi/180   0];
                    [120*pi/180 230*pi/180  90*pi/180];
                    [60*pi/180  60*pi/180   0];
                    [30*pi/180  90*pi/180   120*pi/180]};

% Create basis set:
bas.formalism='sphten-liouv';

% no approximations, i.e complete basis set
bas.approximation='none';

% Use all parts of the coupling tensors to build the spin interaction graph
bas.connectivity='full_tensors';

% Create spin system:

% Create spin system data structure 
spin_system=create(sys,inter);

% Add basis set information
spin_system=basis(spin_system,bas);

% Frame 
spin_system=assume(spin_system,'nmr');

% Liouvillian, states and operators
[H,Q]=hamiltonian(spin_system,'comm');
    H=H+orientation(Q,[alpha beta gamma]);
R=relaxation(spin_system,[alpha beta gamma]);
L=H+1i*R;
kuprov
Posts: 72
Joined: Mon Mar 29, 2021 4:26 pm

Re: Error using alpha

Post by kuprov »

You have not yet defined the three Euler angles that appear in the last few lines. When you request a Hamiltonian, you must specify system orientation - those three alpha, beta, gamma angles in radians.
andrea_simion
Posts: 11
Joined: Sun Aug 15, 2021 9:42 pm

Re: Error using alpha

Post by andrea_simion »

I already specified the euler angles for the isotopes by using inter.zeeman.euler as below, but doesn't work. I have the following error: "Error using alpha. Too many output arguments."

inter.zeeman.euler={[0 0 0];
[0 20*pi/180 10*pi/180];
[90*pi/180 80*pi/180 60*pi/180];
[120*pi/180 130*pi/180 50*pi/180];
[50*pi/180 72*pi/180 0];
[120*pi/180 230*pi/180 90*pi/180];
[60*pi/180 60*pi/180 0];
[30*pi/180 90*pi/180 120*pi/180]};
kuprov
Posts: 72
Joined: Mon Mar 29, 2021 4:26 pm

Re: Error using alpha

Post by kuprov »

Yes, but look at the last three lines of that code. What you had specified is the orientation of interactions in the molecular frame. Now you also need to specify the orientation of the molecule.
andrea_simion
Posts: 11
Joined: Sun Aug 15, 2021 9:42 pm

Re: Error using alpha

Post by andrea_simion »

Thanks, I think I understand now what was wrong. I specified the Euler angles for the interactions in the molecular frame, which in other programs for simulations like spinev or simpson are specified in the CSA file. Now, I need to import a file “rep.100”, where I have 100 different values for alpha, beta, and gamma, and these are used to specify the orientation of the molecule. But the question is: how I can import this .txt file, and say to spinach “Hey Spinach, these are the Euler angles”?
kuprov
Posts: 72
Joined: Mon Mar 29, 2021 4:26 pm

Re: Error using alpha

Post by kuprov »

You don't need to do that; Spinach already has the REPULSION100 grid (see kernel/grids), and also here:

http://spindynamics.org/wiki/index.php? ... wder_grids

If you would like to do the averaging manually, just make a loop, but a better way is to call the powder context (see solid state NMR examples) because that context automatically makes the Hamiltonians, rotates them according to the grid, runs your sequence, and returns the powder average:

http://spindynamics.org/wiki/index.php? ... l_contexts
andrea_simion
Posts: 11
Joined: Sun Aug 15, 2021 9:42 pm

Re: Error using alpha

Post by andrea_simion »

I understand now! Thank you very much for your answers!
I write: parameters.grid='rep_3ang_100pts'; for example, and then I can apply powder context, or if I want to apply Floquet theory, I obtain the Floquet Hamiltonian and then the fid if I call: fid=floquet(spin_system,pulse_sequence,parameters,assumptions);
Now, if I want to apply a specific decoupling sequence to my system (which is not included in the experiments folder of Spinach), e.g tppm, to obtain the fid, I need to call: fid=floquet(spin_system,@tppm,parameters,'nmr');
But first, I must write a code for tppm decoupling including all the details about this pulse sequence. I will try to do that!
kuprov
Posts: 72
Joined: Mon Mar 29, 2021 4:26 pm

Re: Error using alpha

Post by kuprov »

Correct. :) Let me know how you get on.

It is rarely necessary to have three-angle grids, it is likely that a two-angle grid would suffice.
andrea_simion
Posts: 11
Joined: Sun Aug 15, 2021 9:42 pm

Re: Error using alpha

Post by andrea_simion »

It's done :)
I wrote the code for tppm decoupling (see below). I started with the consistency check, then I built the Liouvillian. I defined the initial state, which is Lx on channel 1 (e.g 1H), the detection state, L+ on channel 2 (e.g 13C), and the pulse operators. Then, I started to write the code for the two pulses applied in tppm decoupling, both applied on channel 1 (1H), but with different phases, i.e 0 and 15 degrees. Because in real NMR experiments the power and the duration of rf pulse applied for decoupling are different and need optimizations, I declared these two like parameters (and not as numerical values), i.e parameters.rf_pwr_tppm and parameters.rf_dur_tppm. I split the density operator into two parts, i.e sine and cosine, in each step. I did that also for the fid, but finally I put it together, i.e I wrote the fid like a combination of sine and cosine parts.
There is something that I missed putting in this code? Can you take a look at this and give me suggestions where I need to modify something? There are other optimizations that I need to do?

Code: Select all


% Two pulse phase modulated (tppm) decoupling
%
% Parameters:
%
% parameters.rf_pwr_tppm  - rf power used for tppm decoupling (Hz)
%
% parameters.rf_dur_tppm  - rf pulse duration used for tppm decoupling
%
% parameters.sweep  - sweep widths (Hz)
%
% parameters.npoints  - number of points in fid
%
% parameters.spins - nuclei 
%
% parameters.coil - detection state
%
%    H - Hamiltonian matrix, received from context function
%
%    R - relaxation superoperator, received from context function
%
%    K - kinetics superoperator, received from context function
%
% Outputs:
%
% fid.sin 
% fid.cos
% fid          - sine and cosine components of the State Quadratures and 
%                the total fid (combination of both sine and cosine part)
%
% simi.andrea97@gmail.com
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function fid=tppm(spin_system,parameters,H,R,K)

% Consistency check
grumble(spin_system,parameters,H,R,K);

% Compose Liouvillian
L=H+1i*R+1i*K;

% Get evolution timestep
timestep=1./parameters.sweep;

% Initial state (i.e Lx on spin 1, e.g on 1H)
rho=(state(spin_system,'L+',parameters.spins{1})+state(spin_system,'L-',parameters.spins{1}))/2;

% Detection state (i.e L+ on spin 2, e.g on 13C)
coil=state(spin_system,'L+',parameters.spins{2});

% Pulse operators
% Pulse operator on channel 1 (e.g on 1H)
Lp=operator(spin_system,'L+',parameters.spins{1});
Lx=(Lp+Lp')/2;
Ly=(Lp-Lp')/2i;
% Pulse operator on channel 2 (e.g on 13C)
Fp=operator(spin_system,'L+',parameters.spins{2});
Fx=(Fp+Fp')/2;
Fy=(Fp-Fp')/2i;

% Pulse sequence:

% Step 1: pulse with phase 0 and power parameters.rf_pwr_tppm applied on 
% channel 1 (e.g 1H) for a time parameters.rf_dur_tppm
rho_cos_1=evolution(spin_system,L+2*pi*parameters.rf_pwr_tppm*cosd(0)*Lx,[],rho,parameters.rf_dur_tppm,1,'final');
rho_sin_1=evolution(spin_system,L+2*pi*parameters.rf_pwr_tppm*sind(0)*Ly,[],rho,parameters.rf_dur_tppm,1,'final');

% Step 2: pulse with phase 15 and power parameters.rf_pwr_tppm applied on 
% channel 1 (e.g 1H) for a time parameters.rf_dur_tppm
rho_cos_2=evolution(spin_system,L+2*pi*parameters.rf_pwr_tppm*cosd(15)*Lx,[],rho_cos_1,parameters.rf_dur_tppm,1,'final');
rho_sin_2=evolution(spin_system,L+2*pi*parameters.rf_pwr_tppm*sind(15)*Ly,[],rho_sin_1,parameters.rf_dur_tppm,1,'final');

% Output: total fid
fid.cos=evolution(spin_system,L,parameters.coil,rho_cos_2,timestep(2),parameters.npoints(2)-1,'observable');
fid.sin=evolution(spin_system,L,parameters.coil,rho_sin_2,timestep(2),parameters.npoints(2)-1,'observable');
fid=fid.cos+i*fid.sin;

% Consistency enforcement
    function grumble(spin_system,parameters,H,R,K)
if ~ismember(spin_system.bas.formalism,{'sphten-liouv'})
    error('this function is only available for sphten-liouv formalism.');
end
if (~isnumeric(H))||(~isnumeric(R))||(~isnumeric(K))||(~ismatrix(H))||(~ismatrix(R))||(~ismatrix(K))
    error('H, R and K arguments must be matrices.');
end
if (~all(size(H)==size(R)))||(~all(size(R)==size(K)))
    error('H, R and K matrices must have the same dimension.');
end
if ~isfield(parameters,'sweep')
    error('sweep width should be specified in parameters.sweep variable.');
elseif numel (parameters.sweep)~=2
    error('parameters.sweep array should have exactly two elements.');
end
if ~isfield(parameters,'spins')
    error('working spins should be specified in parameters.spins variable.');
elseif numel (parameters.spins)~=2
    error('parameters.spins cell array should have exactly two elements.');
end
if ~isfield(parameters,'npoints')
    error('number of points should be specified in parameters.npoints variable.');
elseif numel (parameters.npoints)~=2
    error('parameters.npoints array should have exactly two elements.');
end
    end

end

kuprov
Posts: 72
Joined: Mon Mar 29, 2021 4:26 pm

Re: Error using alpha

Post by kuprov »

For one-off pulses, step() is much faster than evolution():

http://spindynamics.org/wiki/index.php?title=Step.m

Also, it is best to keep the grumbler outside the sequence function. At the moment, it is inside.
Post Reply