Calculate angle between vectors

27 views (last 30 days)
Alessandro Ruda
Alessandro Ruda on 4 Mar 2022
Answered: Alessandro Ruda on 9 Mar 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 Comments
Alessandro Ruda
Alessandro Ruda on 6 Mar 2022
Yes, thank you. That is how I should have formulated the question

Sign in to comment.

Answers (2)

Jan
Jan on 4 Mar 2022
Edited: Jan on 4 Mar 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 Comments
Alessandro Ruda
Alessandro Ruda on 6 Mar 2022
Edited: Alessandro Ruda on 6 Mar 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

Sign in to comment.


Alessandro Ruda
Alessandro Ruda on 9 Mar 2022
Anybody?

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by