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

## Trying to simulate DARR

### Re: Trying to simulate DARR

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.

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.