Calculate angle between vectors
28 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Jan
le 4 Mar 2022
Modifié(e) : Jan
le 4 Mar 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.
Réponses (2)
Jan
le 4 Mar 2022
Modifié(e) : Jan
le 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 commentaires
Voir également
Catégories
En savoir plus sur General Applications dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!