Calculate angle between vectors

27 Ansichten (letzte 30 Tage)
Alessandro Ruda
Alessandro Ruda am 4 Mär. 2022
Beantwortet: Alessandro Ruda am 9 Mär. 2022
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
  2 Kommentare
Jan
Jan am 4 Mär. 2022
Bearbeitet: Jan am 4 Mär. 2022
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
Alessandro Ruda am 6 Mär. 2022
Yes, thank you. That is how I should have formulated the question

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Jan
Jan am 4 Mär. 2022
Bearbeitet: Jan am 4 Mär. 2022
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
  5 Kommentare
Alessandro Ruda
Alessandro Ruda am 6 Mär. 2022
Bearbeitet: Alessandro Ruda am 6 Mär. 2022
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
Alessandro Ruda am 6 Mär. 2022
Bearbeitet: Alessandro Ruda am 6 Mär. 2022
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

Melden Sie sich an, um zu kommentieren.


Alessandro Ruda
Alessandro Ruda am 9 Mär. 2022
Anybody?

Kategorien

Mehr zu General Applications finden Sie in Help Center und File Exchange

Produkte


Version

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by