problems with saffron and orientation selection

[JuliaL]

Hello,

I want to do orientation selective HYSCORE simulations using saffron and actually it's two problems I'm encountering...

1)
The first one concerns the definition of atan2 in comparison to eulang. Below the text please find a set of parameters (first two sections) for an electron and a 2H nucleus, with g, A and direction cosines for each. Using both zyz and zxz convention I calculated the Euler angles with acos and atan2 and compared them to the values I got from Mathematica with the same calculations (note: arguments in ArcTan are in reverse order!).
For zyz I also compared them to the results of eulang.
Weird: For only one constellation I had to switch the arguments in Matlab's atan2 to match the control values: Hy_A_psi in zyz.

Is there something in the definitions of atan2 and eulang that can result in this different result? How can I make sure I don't have to check every matrix by hand?


2)
That's the actual saffron problem.
I want to reproduce some older simulations and don't know how EasySpin is internally processing g/A and gpa/Apa. Therefore I want to calculate the orientations via zxz convention on my own and hand over 3x3 g and A matrices to saffron.
But it doesn't not except full g or A matrices. The working example is also below, third section (all together it's an m-file...).
- Giving g and/or A in the form of [x, y, z] and gpa (Apa) as Euler angles, everything is fine.
- Changing g or A to a full 3x3 matrix and keeping gpa/Apa, I get the following error. That's fine, too.

"Error using saffron (line 252)
Full hyperfine matrices are specified in A, so Apa is not allowed!
Error in /Users/.../easyspin-4.5.5/easyspin/saffron.p>saffron (line 76)
Error in JL_saffron_orient (line 85)
saffron(Sys,Exp,Opt)"

- But when I don't use gpa and Apa with the full matrices, I get:

"Error using *
Inner matrix dimensions must agree.
Error in /Users/.../easyspin-4.5.5/easyspin/saffron.p>saffron (line 996)
Error in /Users/.../easyspin-4.5.5/easyspin/saffron.p>saffron (line 76)
Error in JL_saffron_orient (line 85)
saffron(Sys,Exp,Opt)"

I checked the matrix dimensions, took a look on the homepage for the right format. It's 3x3 as indicated.

Please help

Specifications: EasySpin 4.5.5 with Matlab 2012b, IntelMac 10.8.5.

Thanks and kind regards,
Julia (FU Berlin)

Code: Select all

%matlab atan2(y,x) and mathematica arctan[x,y] order reversed!

%Ni, g-values and direction cosines
Ni_g = [2.197 2.139 2.015];
Ni_g_cos = [0.461 0.455 -0.762; -0.204 -0.781 -0.590; -0.864, 0.428, -0.267];


% Euler-angles by hand, zyz-convention, all fine
Ni_g_theta = acos(Ni_g_cos(3,3));
Ni_g_phi = atan2(Ni_g_cos(3,2), Ni_g_cos(3,1));
Ni_g_psi = atan2(Ni_g_cos(2,3), -Ni_g_cos(1,3));

Ni_g_eu = [Ni_g_phi, Ni_g_theta, Ni_g_psi]; 
Ni_g_eulang = eulang(Ni_g_cos);
% compare to Mathematica: 
math_g_eu = [2.6817 1.8411 -0.6589];


% Hydride, hyperfine values and direction cosines
Hy_A = [18.4 -10.8 -18.0]/6.514;
Hy_A_cos = [-0.63 0.54 -0.55; 0.21 0.81 0.55; 0.74 0.23 -0.63];

% Euler-angles by hand, zyz-convention
Hy_A_theta = acos(Hy_A_cos(3,3));
Hy_A_phi = atan2(Hy_A_cos(3,2), Hy_A_cos(3,1));
Hy_A_psi = atan2(-Hy_A_cos(1,3), Hy_A_cos(2,3)); %%%% switched to be identical with eulang %%%%

Hy_A_eu = [Hy_A_phi, Hy_A_theta, Hy_A_psi];
Hy_A_eulang = eulang(Hy_A_cos);
% compare to Mathematica: 
math_A_eu = [0.3013 2.2524 0.7854];


%% zxz-convention, all fine, nothing has to be exchanged

Ni_g_theta_x = acos(Ni_g_cos(3,3));
Ni_g_phi_x = atan2(Ni_g_cos(3,1), -Ni_g_cos(3,2));
Ni_g_psi_x = atan2(Ni_g_cos(1,3), Ni_g_cos(2,3)); 

Ni_g_eu_x = [Ni_g_phi_x, Ni_g_theta_x, Ni_g_psi_x]
% compare to Mathematica: 
math_g_eu_x = [-2.03073 1.8411 -2.22966]

Hy_A_theta_x = acos(Hy_A_cos(3,3));
Hy_A_phi_x = atan2(Hy_A_cos(3,1), -Hy_A_cos(3,2));
Hy_A_psi_x = atan2(Hy_A_cos(1,3), Hy_A_cos(2,3));

Hy_A_eu_x = [Hy_A_phi_x, Hy_A_theta_x, Hy_A_psi_x]
% compare to Mathematica: 
math_A_eu_x = [1.87214 2.2524 -0.785398]

%% Simulation HYSCORE

Exp.Sequence = 'HYSCORE';
Exp.dt = [0.024, 0.024]; %increment
Exp.Flip = [1 1 2 1];
Exp.Inc = [0 1 2 0];
Exp.tau = 0.20; %d_1
Exp.t1 = 0.2;
Exp.t2 = 0.3;
Exp.nPoints = 256;
Exp.Orientations = []; 
Exp.mwFreq = 9.72; %EPR frequency in Ghz
Exp.ExciteWidth = 1000/16;
Exp.Field = 317.0;

Sys.S = 1/2;
Sys.Nucs = '2H';
Sys.Q = [0.03 .1];

Sys.g = Ni_g;
%Sys.g = diag(Ni_g); %%%% does not work %%%%
%Sys.g = [2.197 0 0; 0 2.139 0; 0 0 2.015]; %%%% neither %%%%
Sys.gpa = Ni_g_eu;

Sys.A = Hy_A;
%Sys.A = diag(Hy_A);  %%%% does not work %%%%
%Sys.A = [18.4 0 0; 0 -10.8 0; 0 0 -18.0]/6.511; %%%% neither %%%%
Sys.Apa = Hy_A_eu;

Opt.nKnots = 50;
Opt.ZeroFillFactor = 2;

saffron(Sys,Exp,Opt)

Edit: Similar error occurs on up-to-date stable Debian, Matlab R2011a, Easyspin 4.5.5:
Error using ==> mtimes
Inner matrix dimensions must agree. [...]

[Stefan Stoll]

Hi,
concerning your Euler angle issue, it works fine for me. The second example, where you had to switch atan2 arguments (and the sign), works fine with unswitched arguments. Look here:

Code: Select all

    % Hydride, hyperfine values and direction cosines
    Hy_A = [18.4 -10.8 -18.0]/6.514;
    Hy_A_cos = [-0.63 0.54 -0.55; 0.21 0.81 0.55; 0.74 0.23 -0.63];

    % Euler-angles by hand, zyz-convention
    Hy_A_theta = acos(Hy_A_cos(3,3));
    Hy_A_phi = atan2(Hy_A_cos(3,2), Hy_A_cos(3,1));
    %Hy_A_psi = atan2(-Hy_A_cos(1,3), Hy_A_cos(2,3)); %%%% switched to be identical with eulang %%%%
    Hy_A_psi = atan2(Hy_A_cos(2,3), Hy_A_cos(1,3)); %%%% no switching needed %%%%

    Hy_A_eu = [Hy_A_phi, Hy_A_theta, Hy_A_psi];
    Hy_A_eulang = eulang(Hy_A_cos)
    % compare to Mathematica:
    math_A_eu = [0.3013 2.2524 0.7854]

The output is

Code: Select all

eulang: Rotation matrix is not orthogonal within 1e-6.
Hy_A_eulang =
    0.3040    2.2492    0.7852
math_A_eu =
    0.3013    2.2524    0.7854

The difference between eulang() and your manual calculation is that eulang() actually uses a least-squares fitting approach that minimizes the error between your input rotation matrix and the one calculated using erot() and the angles.

You are giving a rotation matrix with only 2 significant figures in each element. This is fairly inaccurate. Depending on the way you calculate the angles, you will get results that differ by 1 degree or so. eulang() avoids these problems by the least-squares fitting procedure.

[Stefan Stoll]

Re Hyscore problem: This might be a bug. Can you post a minimal example that shows the bug in a separate post? Thanks

[JuliaL]

Hallo Stefan,

thanks for the answer!
I don't know why, but for me the switch seems essential. Without:

Code: Select all

Hy_A_eu = 0.3013    2.2523    2.3562
eulang: Rotation matrix is not orthogonal within 1e-6.
Hy_A_eulang = 0.3040    2.2492    0.7852
math_A_eu = 0.3013    2.2524    0.7854

[Stefan Stoll]

With ES 4.5.5, with

Code: Select all

R = [-0.63 0.54 -0.55; 0.21 0.81 0.55; 0.74 0.23 -0.63];
abc = eulang(R)
R1 = erot(abc);
norm(R-R1)

I get the following output:

Code: Select all

eulang: Rotation matrix is not orthogonal within 1e-6.
abc =
    0.3040    2.2492    0.7852
ans =
    0.0064

What do you get?

[JuliaL]

I get the same.

Concerning your first answer: Hy_A_eu = [Hy_A_phi, Hy_A_theta, Hy_A_psi]; is commented, you just show the results by eulang. Which value do you have for Hy_A_psi?

[Stefan Stoll]

The correct way to calculate the third Euler angle from the rotation matrix is

Code: Select all

Hy_A_psi = atan2(Hy_A_cos(2,3), -Hy_A_cos(1,3))

That gives 0.7854.

[JuliaL]

Found a missing sign in my second to last post but that wasn't the problem in the first place.
I don't know why it works know, I didn't change anything. But MatLab messes with my stuff from time to time, anyways.

Thanks again for the effort, thanks for the explanation of how eulang works! Next time I'll keep on repeating until it works!

(I noticed, that for this special case atan2(-Hy_A_cos(1,3), Hy_A_cos(2,3)) and atan2(Hy_A_cos(2,3), -Hy_A_cos(1,3)) give the same result.)