Trying to simulate DARR

Topics related to Spinach package
Post Reply
rgpg05
Posts: 1
Joined: Thu Nov 03, 2022 6:55 pm

Trying to simulate DARR

Post by rgpg05 »

Hello Prof. Illya

I am trying to simulate DARR polarization transfer between two carbons in the presence of 4 protons (to start with). I have written the code. I tested this for some very simple cases (like rotational resonance by removing protons and setting rf to zero) and it seems to work fine. On including more than two protons, simulations slow down considerably (Note : no basis approx./ cut off tolerance used at this point). Is there something, I am doing terribly wrong. I would appreciate if you could suggest some ways to improve on its performance.

I also had a query regarding the timestep over which Hamiltonian is considered constant. In SIMPSON, it is 'maxdt'. What is the analogue in SPINACH.

Below is the code.



*******************************************MAIN CODE*********************************
function darr_pdsd()

% System specifications
[sys,inter]=s2spinach('spinsys.in');
sys.magnet=14.08;


% Basis set specifications
bas.formalism='sphten-liouv';
bas.approximation='none';
bas.connectivity='full_tensors';

% Spinach housekeeping
spin_system=create(sys,inter);
spin_system=basis(spin_system,bas);

% MAS parameters
parameters.rate = 11000.00;
parameters.axis = [1 1 1 ];
parameters.max_rank=10;
parameters.grid='rep_2ang_100pts_sph';
% density, detection and active channels
% factor of 4 included to normalize magnetization to unity
parameters.spins={'13C','1H'};
parameters.rho0=state(spin_system,{'Lz'},{1});
parameters.coil=4.0*state(spin_system,{'Lz'},{2});
%additional parameters
parameters.verbose=1;
parameters.np=1000;
parameters.dw=1e-5;
parameters.rf=parameters.rate;
% Simulation
fid=singlerot(spin_system,@darr,parameters,'nmr');
% Plot the answer
x_axis=linspace(0,parameters.dw*parameters.np,parameters.np+1);
figure(); plot(x_axis,real(fid)); kgrid;
ylabel('build up on second carbon');
xlabel('mixing time, seconds');
end
******************************************PULSE SEQUENCE*********************************
function fid=darr(spin_system,parameters,H,R,K)
% Check consistency
grumble(parameters,H,R,K);
% Compose the Liouvillian
L=H+1i*R+1i*K;
% Irradiation on 1H channel
Hp=operator(spin_system,'L+',parameters.spins{2});
Hp=kron(speye(parameters.spc_dim),Hp);
Hx=(Hp+Hp')/2;Hy=(Hp-Hp')/2i;
%add to the Liouvillian
L=L+2*pi*parameters.rf*Hx;
%acquire
fid=evolution(spin_system,L,parameters.coil,parameters.rho0,...
parameters.dw,parameters.np,'observable');

end
% Consistency enforcement
function grumble(parameters,H,R,K)
if (~isnumeric(H))||(~isnumeric(R))||(~isnumeric(K))||...
(~ismatrix(H))||(~ismatrix(R))||(~ismatrix(K))
error('H, R and K arguments must be matrices.');
end
if ~isfield(parameters,'np')
error('number of time steps should be specified in parameters.np variable.');
end
if ~isfield(parameters,'rho0')
error('initial state must be specified in parameters.rho0 variable.');
end
if ~isfield(parameters,'coil')
error('detection state must be specified in parameters.coil variable.');
end
end
******************************************************************************************
thanks for your time

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

Re: Trying to simulate DARR

Post by kuprov »

Dear Rajat,

of course it would slow down as you add more spins - without approximations, the complexity of Liouville space simulations of spin-1/2 systems grows very rapidly. Matrix dimensions go as 4^n, and then matrix-vector multiplication complexity scales quadratically with the dimension, so we get 8^n complexity scaling. I am afraid ten spins is about the limit. Beyond that, you need restricted state spaces.

Spinach does not have a built-in time step variable; the time step is set by the pulse sequence (meaning, by whoever wrote that) or simply by the user. For evolution without detection, it's arbitrary (so long as you evolve for the required total time). For evolution with detection, it is determined by the Nyquist-Shannon condition: it is necessary to have two points per period of the fastest frequency.

Your code is pretty much optimal, make sure the calculation is converged with respect to the powder grid and the rotor rank, rest is fine.

To use restricted state spaces, follow the guidance here (https://spindynamics.org/wiki/index.php ... cification). With more than ten spins, that is more or less unavoidable.

Best wishes,
Ilya.
Post Reply