Calculate angle between vectors
27 views (last 30 days)
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.
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)];