Two spin systems in a custom function

[thanasis]

I have written a script that fits pairs of spins that interact through dipolar interaction using the model in T.D. Smith, J.R. Pilbrow, Coord. Chem. Rev. 1974, 13, 173-278 (http://dx.doi.org/10.1016/S0010-8545(00)80255-6), equation 14.

Since the interaction energy is a function of an external variable, the distance r, I introduced that into my system structure, along with the ξ and η angles. Then I created a system of two S=1/2 spins (Sys.S=[1/2 1/2]) and used a custom function to calculate the interaction matrix as is described in the paper. This function also fixes their gx/y/z and their gStrainx/y/z parameters to the same values during fitting, since they are supposed to be similar. It took some time to type in the matrix elements, but it works very well and the results make sense.

Now, I want to introduce a second system, comprising a single spin. The reason I need a separate system is that the dipole-inetracting pair and the isolated spin have different relaxation times, i.e. different intensities, and I need to adjust this through their weights in Easyspin. So I made Sys1.S = [1/2 1/2] and Sys2.S = 1/2, but I have trouble introducing both Sys1 and Sys2 into my custom function.

When I try:

Code: Select all

function [x,y] = custom({Sys1,Sys2},Exp,Opt)
fullSys = {Sys1,Sys2};
etc...

It seems that the custom function is not ready to accept two systems:

Error using nargout
Error: File: custom.m Line: 1 Column: 38
Unbalanced or unexpected parenthesis or bracket.
Error in esfit (line 151)
Error in S12_S12_S12_dipolar (line 77)
esfit('custom',spc,{Sys1,Sys2},{Vary1,Vary2},Exp,[],FitOpt);

Is my approach valid for Easyspin? Should I use another approach, e.g. calculate the dipole-pair matrix in my custom function first and then include the second system in the main script?

Thanks for the help!

[Matt Krzyaniak]

Your approach should be valid, esfit loops over the number of elements in the sys cell you provide and passes them each to the function, so likely what you need is an if/else statement in your custom function.


for instance

Code: Select all

function [x,y]=custom(sys,exp,opt)
if length(sys.S)==1
[x,y]=pepper(sys,exp,opt)
...
else
twospin interaction stuff
...
end

[thanasis]

Thanks Matt.

In this approach, the "sys" that you propose in the function call (function [x,y]=custom(sys,exp,opt)) can it be something generic, or something that has already been defined in my main script?

Shouldn't I make specific references to Sys1 (dimer system) and Sys2 (monomer system) which I have defined in my main script?

E.g. like

Code: Select all

function [x,y]=custom(sys,exp,opt)

if length(sys.S)==2
[x,y]=pepper(Sys1,exp,opt)
...interaction...

if length(sys.S)==1
[x,y]=pepper(Sys2,exp,opt)
...monomer...

end

Is this what it should look like?

[thanasis]

I have also tried some alternatives, having defined my Sys1 (pair) and Sys2 (monomer):

Code: Select all

Sys1.S = [1/2 1/2]
...

Sys2.S = 1/2
...

Vary1...
Vary2...

Sys = {Sys1,Sys2};
Vary = {Vary1,Vary2};

esfit('custom',spc,Sys,Vary,Exp,[],FitOpt);
...

and in the custom function:

Code: Select all

function [x,y] = custom(Sys,Exp,Opt)

if length(Sys)==2 % Two spins
fullSys = Sys{1,1}; % First array in Sys is Sys1
fullSys.g = [Sys{1,1}.g];
fullSys.ee = dipolar... etc
[x,y] = pepper(fullSys,Exp,Opt);

else
else
fullSys = Sys{1,2};
[x,y] = pepper(fullSys,Exp,Opt);
end

return

I manage to get the fitting GUI, but when I launch I get the error:

Cell contents reference from a non-cell array object.
Error in custom (line 23)
fullSys = Sys{1,2};

However, I don't get that error for fullSys = Sys{1,1}

If I change the curly to normal brackets: fullSys = Sys(1,2);

Code: Select all

Index exceeds matrix dimensions.
Error in custom (line 23)
fullSys = Sys(1,2);

It seems I can't successfully invoke the components Sys1 and Sys2 of the composite Sys structure from my custom function.

[Matt Krzyaniak]

First to clarify how esfit works, if you provide it a multi-structure cell input, sort of like how you're trying to set it up now, esfit will loop over the individual components and pass them to the function so in terms of code it looks something like:

Code: Select all

fit = 0;
for i=1:length(sys)
[~,y] = custom(sys{i},exp,opt);
fit = fit+y;
end

So if you write your custom function to try and separate out the components in the input cell you're bound to run into problems. I was suggesting that you write your custom function in a general way and use the if/elseif/else to determine how you treat the specific input the function receives. So you where almost there in your first reply. it should have been:

Code: Select all

function [x,y]=custom(sys,exp,opt)

if length(sys.S)==2
[x,y]=pepper(sys,exp,opt)
...interaction...

elseif length(sys.S)==1
[x,y]=pepper(sys,exp,opt)
...monomer...

end

You have to remember, inside a function you can only access variables that you've passed to the function, so only sys, exp, and opt are available.

[thanasis]

OK, I think I get it a little better about esfit.

What I haven't understood is where to define the individual structures, (Sys1.S, Sys1.g, Sys1.weight,..., and Sys2.S,...). Should I do that in my main function and then pass them on to the custom one, or should I define them in the custom function altogether?

If I define them in my main function, I currently understand that the correct syntax to pass them on is:

Code: Select all

esfit('custom',spc,{Sys1,Sys2},{Vary1,Vary2},Exp,[],FitOpt);

But on the custom function side, what is the "Sys" that I will invoke? What is the correct syntax to pass them on and then invoke them?

Or have I entirely misunderstood?

Thanks for the help!

[Matt Krzyaniak]

For a more general custom function define your system in the main script. I'd set things up like this:

Code: Select all

% parameter testing script
clear sys1 sys2 exp opt
sys1.S = 1/2;
sys1.g = 2;
....
sys1.weight = 1;

sys2.S = [1/2 1/2];
sys2.g = [2 2];
...
sys2.weight = 1;

exp.Range = [300 400];
exp.nPoints = 1024;
exp.mwFreq = 9.5;

opt  = []
[x,y1] = custom(sys1,exp,opt);
[~,y2] = custom(sys2,exp,opt);
plot(x,y1,x,y2,x,y1+y2)
%% fitting

vary1.g = 0.1;
vary2.g = [0.1 0.1];
esfit('custom',spc,{sys1,sys2},{vary1,vary2},Exp,[],FitOpt);

with a custom function:

Code: Select all

function [x,y]=custom(sys,exp,opt)

if length(sys.S)==2
[x,y]=pepper(sys,exp,opt)
...interaction...

elseif length(sys.S)==1
[x,y]=pepper(sys,exp,opt)
...monomer...

end

[thanasis]

Yes, that's it!

I got it to work as you suggested and it seems correct.

Thanks a lot!

[thanasis]

I am resurrecting this thread as it seems esfit has changed the way it works in the meantime (I am using ES 6.0.0-dev.10).

I have defined Sys1 and Sys2, and have a custom function custom_2sys, which I invoke for esfit as follows:

Code: Select all

esfit('custom_2sys',spc,{Sys1,Sys2},{Vary1,Vary2},Exp,[],FitOpt);

Then, in that custom function I try to decide which system it is currently calculating, by using Matt's previous method:

Code: Select all

if Sys.r12 == 3.486 % First complex
    fullSys.g = [Sys1.gxy Sys1.gxy Sys1.gz; Sys1.gxy Sys1.gxy Sys1.gz; Sys1.gxy Sys1.gxy Sys1.gz];
elseif Sys.r12 == 3.410 % Second complex
    fullSys.g = [Sys2.gxy Sys2.gxy Sys2.gz; Sys2.gxy Sys2.gxy Sys2.gz; Sys2.gxy Sys2.gxy Sys2.gz];
end

or

Code: Select all

if Sys.r12 == 3.486 % First complex
    fullSys.g = [Sys{1,1}.gxy Sys{1,1}.gxy Sys{1,1}.gz; Sys{1,1}.gxy Sys{1,1}.gxy Sys{1,1}.gz; Sys{1,1}.gxy Sys{1,1}.gxy Sys{1,1}.gz];
elseif Sys.r12 == 3.410 % Second complex
    fullSys.g = [Sys{1,2}.gxy Sys{1,2}.gxy Sys{1,2}.gz; Sys{1,2}.gxy Sys{1,2}.gxy Sys{1,2}.gz; Sys{1,2}.gxy Sys{1,2}.gxy Sys{1,2}.gz];
end

Both get me an error:

Dot indexing is not supported for variables of this type.

It seems like esfit now imports both system structures simultaneously as cells, so Matt's trick doesn't work.
I can write my if statement like if Sys{1,1}.r12 == 3.486, to get around this error, but that would beat the whole purpose, wouldn't it?
So how can I know which system esfit is currently calculating?
The other option would be to rewrite my custom function to explicitly calculate each system. But as it is already quite big, it would be an issue to rewrite all the code.

[Stefan Stoll]

Please update to the latest 6.0.0-dev version, since functions such as esfit are continuously improving.

esfit now calls the custom simulation function with a cell array of all spin systems, i.e. custom_2sys({Sys1,Sys2},.... So all you need to do is to write a loop inside the simulation function.