# Calculate angle between vectors

27 Ansichten (letzte 30 Tage)
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):
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 KommentareKeine anzeigenKeine ausblenden
Jan am 4 Mär. 2022
Bearbeitet: Jan am 4 Mär. 2022
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 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 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 Kommentare3 ältere Kommentare anzeigen3 ältere Kommentare ausblenden
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 am 6 Mär. 2022
Bearbeitet: Alessandro Ruda am 6 Mär. 2022
Ultimately, this is the code:
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.');
H3H1p_angle_rad = atan2(vecnorm(myCross2(H3H1p,Z), 2, 2), H3H1p * Z.');
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 am 9 Mär. 2022
Anybody?
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

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

R2020b

### Community Treasure Hunt

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

Start Hunting!

Translated by