Calculate angle between vectors

Dear MatLab comunity,

I have a vectorial problem. I am writing again the question because it wasn't clearly explained, so I'll try to be as precise as possible.

I have a molecular dynamics trajectory of a molecule with 1000 frames.

For this molecule I calculated the three major axis of inertia (see picture)

The three major axes of inertia intersect at the center of mass of coordinates:

COM = [3.9721615314483643 -8.11227798461914 -0.34564414620399475]

and each one has coordinates, respectively:

X_ax = [0.43753358721733093 0.719929575920105 -0.5387632250785828] 
Y_ax = [0.5724487900733948 0.2390475869178772 0.7843204736709595] 
Z_ax = [0.6934455633163452 -0.6515809297561646 -0.3075314462184906]

Now, I have chosen three protons in the molecule: H1p, H5p and H3.

Of these I extracted the xyz coordinates in each frame and stored them (see file attached, the first column is the frame number):

L_H1p = load('H1p.dat');
L_H5p = load('H5p.dat');
L_H3 = load('H3.dat');

H1p = [L_H1p(:,2), L_H1p(:,3), L_H1p(:,4)];
H5p = [L_H5p(:,2), L_H5p(:,3), L_H5p(:,4)];
H3 = [L_H3(:,2), L_H3(:,3), L_H3(:,4)];

so having an array of coordinates for each atom.

Now, I want to define the vectors that goes from H1p to H5p and the vector that goes from H1p to H3 and I did it as following:

H5pH1p = [(H5p(:,1)-H1p(:,1)), (H5p(:,2)-H1p(:,2)), (H5p(:,3)-H1p(:,3))];
H3H1p = [(H3(:,1)-H1p(:,1)), (H3(:,2)-H1p(:,2)), (H3(:,3)-H1p(:,3))];

Now I have the array of 1000 vectors for H5pH1p and 1000 vectors for H3H1p.

I want to study how these vectors move with reference to the axis of inertia.

Let's say that I want to calculate the angle between the array of vectors H5pH1p and the Z component of the inertia axes.

I am not sure how to approach to this and hopefully I stated clearly the problem.

Best,

Alex

[Jan]

As mentioned in another thread already: You can simplify

H5pH1p = [(H5p(:,1)-H1p(:,1)), (H5p(:,2)-H1p(:,2)), (H5p(:,3)-H1p(:,3))];

to

H5pH1p = H5p - H1p;

or:

H1p = [L_H1p(:,2), L_H1p(:,3), L_H1p(:,4)];

to

H1p = L_H1p(:,2:4);

This reduces the cance of typos and is much faster.

The actual question is: "I want to calculate the angle between the array of vectors H5pH1p and the Z component of the inertia axes" . Then all you need to know to create an answer is the dimension of the vectors. So "X = rand(1000, 3), Z = rand(1, 3)" is sufficient and much shorter than your question.

[Alessandro Ruda]

Yes, thank you. That is how I should have formulated the question

Answer 1

[Jan]

See: https://www.mathworks.com/matlabcentral/answers/16243-angle-between-two-vectors-in-3d

The actual numerically stable formula is: atan2(norm(cross(X, Z)), dot(X, Z));

Unfortunately Matlab's cross and dot do not work with implicit expanding yet.

X = rand(1000, 3);
Z = rand(1, 3);
A = atan2(vecnorm(myCross(X, Z), 2, 2), X * Z.');

function Z = myCross(X, Y)
% INPUT:   X, Y: [M x 3] or [1 x 3] arrays.
% OUTPUT:  Z: Cross product.
Z = [ X(:,2) .* Y(:,3) - X(:,3) .* Y(:,2), ...
      X(:,3) .* Y(:,1) - X(:,1) .* Y(:,3), ...
      X(:,1) .* Y(:,2) - X(:,2) .* Y(:,1)];
end

[Alessandro Ruda]

Sorry, I don't quite get what is happening in the code! I've seen all the other threads and the formulation with arctan is the most stable one, okay.

Then vecnorm returns the magnitude of the cross products between the vectors X and Z and scalar product between the same, before the arctan is calculated.

As you say "Matlab's cross and dot do not work with implicit expanding yet" which i didn't quite get but i guess it means that doesn't handle products of matrixes with different lenght or something so we need to defined the cross product explicitly.

Then you define the cross product function where I guess X is H5pH1p and Y is the Z_axis [1, 3] vector.

But I don't get why if Z is already defined, you define a function Z afterward? sorry, it might be stupid, hope you have patience with this, I am just trying to understand.

I set it this way:

A = atan2(vecnorm(myCross(H5pH1p,Z_ax), 2, 2), H5pH1p * Z_ax.');


function Z = myCross(H5pH1p, Z_ax)
% INPUT:   X, Y: [M x 3] or [1 x 3] arrays.
% OUTPUT:  Z: Cross product.
Z = [ H5pH1p(:,2) .* Z_ax(:,3) - H5pH1p(:,3) .* Z_ax(:,2), ...
      H5pH1p(:,3) .* Z_ax(:,1) - H5pH1p(:,1) .* Z_ax(:,3), ...
      H5pH1p(:,1) .* Z_ax(:,2) - H5pH1p(:,2) .* Z_ax(:,1)];
end

First I am not sure if I have to use the vector Z = Z_ax - COM instead of Z_ax directly.

Second, Is it not that the vectors of the first matrix (H5pH1p) have to be translated to the same origin of the the Z_ax vector for the angle to be calculated and the cross product to be obtained, or am I just confused about it? Because my reference point are the three inertia axes with origins in COM.

Anyway, with the code as it is I get a 1000x1 matrix of angle (in radians I guess). Is this the array of angles between each vector in the array H5pH1p and the Z axis ?

Thank you for the help!

[Jan]

The actual algorithm is:

Angle = atan2(norm(cross(X, Z)), dot(X, Z))

For the implementation in Matlab we have to consider, that norm() calculates the matrix norm, when the argument is a matrix. Therefore vecnorm is the wanted function, which calculates the norm over vectors along the specified direction.

Many Matlab function are working with "implicit expanding". This means, that if the operands are [M x N] and [1 x N], the vector is expanded to the same size automatically. This works e.g. for +, -, .*, ./ and so on, but not for the functions cross() and dot().

"Then you define the cross product function where I guess X is H5pH1p and Y is the Z_axis [1, 3] vector." - almost. It is very useful to implement a function not only for a specific input, but generally. If you call e.g. sin(x), it does not matter, if the argument is called "x", "y", "hurliturli_Kafonz" or if it is constand like 19. Using the names "H5pH1p" and "Z_ax" as inputs is not useful, because "X" and "Y" are short, clear and general.

The function is defined with its inputs and outputs by the line

function Z = myCross(X, Y)

Now the names of the variables in the caller do not matter in any way, because the function sees their values only. The "Z" in the caller has no relation to the "Z" used inside the function. Both are variables and each function has its own set of variables completely independent from other functions.

You can use the function "myCross()" exactly how I have posted it. There is no need to rename the variables.

"First I am not sure if I have to use the vector Z = Z_ax - COM instead of Z_ax directly." - I do not know this either. It depends on what you want to calculate.

"Is it not that the vectors of the first matrix have to be translated to the same origin of the the Z_ax vector for the angle to be calculated and the cross product to be obtained" - as said before, a vector is a set of e.g. 3 elements. You can interprete it as coordinates of a point in 3D or a direction. If you want to calculate angles between vectors, there is no need to move them around in space, if the vectors are directions. If they are coordinates, and "angle" between them is not defined mathematically. So it really depends on the mathematical nature of the problem, but I did not understand yet what you want to calculate.

"Is this the array of angles between each vector in the array H5pH1p and the Z axis" - exactly.

[Alessandro Ruda]

Thank you for the quick reply and the exhaustive explanation!

I am still not sure I am calculating what I have in mind, of course due to my incapability of defining the problem correctly.

All the coordinates I obtained are calculated from an origin point that lies far from the molecule and its center of mass, and that is the general reference system.

What would actualy be conveninent for me is to have the reference system centered where the center of mass of the molecules is i.e. with O(0,0,0) corresponding to the COM = [3.9721615314483643 -8.11227798461914 -0.34564414620399475].

But since that is not the case, my guess is that I have to take into consideration that the vectors H5p - H1p and Z_ax - COM do not share the same origin for the way I am seeing them.

And I guess this is what I didn't make clear.

My reference vector is Z = Z_ax - COM where COM is my origin coordinates, ||Z_ax - COM|| is the magnitude of the vector and Z_ax are the coordinates of the vector that give the direction to it with respect to COM tough, not the real origin O(0,0,0).

Therefore my reference vector is Z = Z_ax - COM, and COM is the origin of my reference system.

The vector H5pH1p is simply the distance between two protons in the molecule with ||H5p-H1p|| being its magnitude.

Now, to calculate the angle between this vector and the Z vector I think I should first translate H5pH1p to the COM origin and only then apply the formula you provided me to calculate the angle between the two.

Otherwise my view is that I am calculating the angle between the vectors H5pH1p - O(0,0,0) and Z - O(0,0,0).

Am I interpreting this correctly?

I hope I am gettin more clear with the explanation, I am still a novice with matlab and vector calculus is something I scratched the surface of at the beginnig of my studies many years ago.

[Alessandro Ruda]

So, I guess what I have to do would be to translate every vector with the vector COM as follow:

A = H1p - COM;
B = H5p - COM;
C = H3 - COM;
Z = Z_ax;

then do the calculation with H5pH1p defined as:

H5pH1p = B - A;

as you told me:

A = atan2(vecnorm(myCross(H5pH1p,Z), 2, 2), H5pH1p * Z.');
function Z = myCross(H5pH1p, Z)
% INPUT:   X, Y: [M x 3] or [1 x 3] arrays.
% OUTPUT:  Z: Cross product.
Z = [ H5pH1p(:,2) .* Z(:,3) - H5pH1p(:,3) .* Z(:,2), ...
      H5pH1p(:,3) .* Z(:,1) - H5pH1p(:,1) .* Z(:,3), ...
      H5pH1p(:,1) .* Z(:,2) - H5pH1p(:,2) .* Z(:,1)];
end

Or am I still wrongly seing it?

/alex

[Alessandro Ruda]

Ultimately, this is the code:

L_H1p = load('H1p.dat');
L_H5p = load('H5p.dat');
L_H3 = load('H3.dat');

H1p = L_H1p(:,2:4);
H5p = L_H5p(:,2:4);
H3 = L_H3(:,2:4);

COM = [3.9721615314483643 -8.11227798461914 -0.34564414620399475]  %center of mass coordinates
X = [0.43753358721733093 0.719929575920105 -0.5387632250785828] 
Y = [0.5724487900733948 0.2390475869178772 0.7843204736709595] 
Z = [0.6934455633163452 -0.6515809297561646 -0.3075314462184906]

% Translate the coordinates to the new origin

H1p_o = H1p - COM;
H5p_o = H5p - COM;
H3_o = H3 - COM;

% Generate vectors

H5pH1p = H5p_o - H1p_o;
H3H1p = H3_o - H1p_o;

% Calculate angles with respect to the Z-axis of inertia

H5pH1p_angle_rad = atan2(vecnorm(myCross(H5pH1p,Z), 2, 2), H5pH1p * Z.');
H5pH1p_angle_deg = H5pH1p_angle_rad*180/pi;

H3H1p_angle_rad = atan2(vecnorm(myCross2(H3H1p,Z), 2, 2), H3H1p * Z.');
H3H1p_angle_deg = H3H1p_angle_rad*180/pi;

function Z = myCross(H5pH1p, Z)
Z = [ H5pH1p(:,2) .* Z(:,3) - H5pH1p(:,3) .* Z(:,2), ...
      H5pH1p(:,3) .* Z(:,1) - H5pH1p(:,1) .* Z(:,3), ...
      H5pH1p(:,1) .* Z(:,2) - H5pH1p(:,2) .* Z(:,1)];
end

function Z = myCross2(H3H1p, Z)
Z = [ H3H1p(:,2) .* Z(:,3) - H3H1p(:,3) .* Z(:,2), ...
      H3H1p(:,3) .* Z(:,1) - H3H1p(:,1) .* Z(:,3), ...
      H3H1p(:,1) .* Z(:,2) - H3H1p(:,2) .* Z(:,1)];
end

Answer 2

[Alessandro Ruda]

Anybody?