Output order for resfreqs_matrix and resfreqs_perturb

[katarkon]

Surprisingly, the order of line positions calculated by resfreqs_matrix() and resfreqs_perturb() is different. This makes the functions non-interchangeable. It should be better that both functions return equal results.
Here the script illustrating the problem.

Code: Select all

Sys3.g = 2.00552;
Sys3.A = mt2mhz([7.82, 1.89]/10); 
Sys3.Nucs = '14N,1H';
Param.Field=344;
Op.PerturbOrder=2;
[pos,int]=resfreqs_matrix(Sys3,Param,Op);
disp(pos');
[pos,int]=resfreqs_perturb(Sys3,Param,Op);
disp(pos');

It returns the follow output:

Code: Select all

>> resfreq_test
    9.6806    9.6587    9.6753    9.6367    9.6534    9.6315

    9.6315    9.6367    9.6534    9.6587    9.6753    9.6806

The line positions are the same, but their order is different.

[Stefan Stoll]

This is indeed an inconsistency. We will open an issue and investigate it. Thanks for reporting!

[katarkon]

I see the problem is still unsolved in last version (5.2.11). Is any progress in the solving?
I found that the problem has unpleasant effect on the order of output for pepper if Opt.Output='separate' is defined. It is absolutely different for matrix, perturbation or hybrid treatment. Moreover, for last method the order changes absolutely unpredictable with including of some nuclei in Opt.HybridCoreNuclei.
I think the problem is the order of transition selection inside the function resfields. Here the script illustrating the problem

Code: Select all

Sys.g = 2;
Sys.Nucs = '14N,14N';
Sys.lwpp = [0,0.02];
Sys.A = [15; 15]; % MHz
Exp.Range=[340 348];
Exp.mwFreq = 9.65;
[Pos,Amp,Wid,Trans,Grad] = resfields(Sys,Exp);
disp(Trans);
[Pos,Amp,Wid,Trans,Grad] = resfields_perturb(Sys,Exp);
disp(Trans);

The attempt to use the list of transitions was failed.

Code: Select all

Sys.g = 2;
Sys.Nucs = '14N,14N';
Sys.lwpp = [0,0.02];
Sys.A = [15; 15]; % MHz
Exp.Range=[340 348];
Exp.mwFreq = 9.65;
[Pos,Amp,Wid,Trans,Grad] = resfields_perturb(Sys,Exp);
disp(Trans);
Opt.Transitions=Trans;
[Pos,Amp,Wid,Trans,Grad] = resfields(Sys,Exp,Opt);
disp(Trans);

This is completely unexpected. At last, if Opt.Transitions='all' is used the results are more strange.

Code: Select all

Sys.g = 2;
Sys.Nucs = '14N,14N';
Sys.lwpp = [0,0.02];
Sys.A = [15; 15]; % MHz
Exp.Range=[340 348];
Exp.mwFreq = 9.65;
Opt.Transitions='all';
[Pos,Amp,Wid,Trans,Grad] = resfields_perturb(Sys,Exp,Opt);
disp(Trans);
[Pos,Amp,Wid,Trans,Grad] = resfields(Sys,Exp,Opt);
disp(Trans);
[Pos,Amp,Wid,Trans,Grad] = resfields_perturb(Sys,Exp);
disp(Trans);
[Pos,Amp,Wid,Trans,Grad] = resfields(Sys,Exp);
disp(Trans);

UPD. Opt.Transitions='all' affects on the output order of pepper function too. The order of the spectra becomes less unpredictable for matrix and hybrid treatments but just differs from perturbation method.

[Stefan Stoll]

This issue is still unresolved. There is a bug in the level indices in the transition output of resfreqs_perturb and resfields_perturb which needs to be fixed first before the sorting can be addressed.

Concerning the sorting, we will likely not be able to guarantee ordering that is identical between all functions and all options settings, but we will try.

[katarkon]

I suggest that sorting of output of the function pepper according to resonance fields values (e.g. in order of increasing the field) should be enough for temporary solution. A special key, such as Opt.Output = 'sorted', can be used to not violate the compatibility.
Fixing the bug with level indexing is better, of course.

[katarkon]

The rough solution for the problem is ordering of the resfields output according to resfields_perturb. Although forbidden transitions could make a certain problem.

Code: Select all

[Pos] = resfields_perturb(Sys,Exp,Opt);
[Pos1] = resfields(Sys,Exp,Opt);
Ii=[];

for i=1:numel(Pos);
    Pos2=abs(Pos1-Pos(i));
    while true
        [C,I]=min(Pos2);
        if sum(Ii==I)<1;
            Ii(i)=I;
            break
        end
        Pos2(I)=max(Pos2);
    end

end
Pos1=Pos1(Ii);

[katarkon]

I found that the QZ algorithm for the function eig does not violate the order of the eigenvalues. This can be a step towards solving the problem. The alternative form eig(H,eye(size(H)),'qz') instead eig(H) may be used for testing.

[Stefan Stoll]

resfields and other functions that are based on spin Hamiltonian matrix diagonalization are designed to handle any spin system, including ones with very strong state mixing. The only ordering method of energy levels that works across all spin systems in a consistent fashion is to order them by increasing energy, which is done by eig. There is no predefined general intrinsic order of energy levels, except when the states are almost pure Zeeman states. However, these represent only a small subclass of possible spin systems.

[katarkon]

Thanks for explanation. The ordering of the levels according to the energy should avoid their crossing, I suppose.