How to calculate torsion angle?

Given the x-, y-, and z-coordinates of four points in 3-dimensional (3D) space, how to calculate torsion angle? Overall, this is a well-solved problem in structural biology, and one can find detailed description of the algorithm in text books and on-line documents. The algorithm for calculating torsion angle is implementated in virtually every software package in structural chemistry or biology.

Basic as it is, however, in my experience, it is very important to have a detailed appreciation of how the method works in order to really get into the 3D world. Here is a worked example using Octave/Matlab of my simplified, geometry-based implementation of how to calculate torsion angle, including how to determine its sign. No theory or (complicated) mathematical formula, just a step-by-step illustration of how I solve this problem.

1. Coordinates of four points in variable abcd:
   

   abcd = [ 21.350  31.325  22.681
   22.409  31.286  21.483
   22.840  29.751  21.498
   23.543  29.175  22.594 ];

2. Two auxiliary functions: norm_vec() to normalize a vector; get_orth_norm_vec() to get the orthogonal component (normalized) of a vector with reference to another vector, which should have been normalized.
   

   function ovec = norm_vec(vec)
     ovec = vec / norm(vec);
   endfunction
   
   function ovec = get_orth_norm_vec(vec, vref)
     temp = vec - vref * dot(vec, vref);
     ovec = norm_vec(temp);
   endfunction

3. Get three vectors: b_c is the normalized vector b→c; b_a_orth is the orthogonal component (normalized) of vector b→a with reference to b→c; c_d_orth is similarly defined, as the orthogonal component (normalized) of vector c→d with reference to b→c.
   

   b_c = norm_vec(abcd(3, :) - abcd(2, :))
     % [0.2703158  -0.9627257   0.0094077]
   b_a_orth = get_orth_norm_vec(abcd(1, :) - abcd(2, :), b_c)
     % [-0.62126  -0.16696   0.76561]
   c_d_orth = get_orth_norm_vec(abcd(4, :) - abcd(3, :), b_c)
     % [0.41330   0.12486   0.90199]

4. Now the torsion angle is the angle between the two vectors, b_a_orth and c_d_orth, and can be easily calculated by their dot product. The sign of the torsion angle is determined by the relative orientation of the cross product of the same two vectors with reference to the middle vector b→c. Here they are in opposite direction, thus the torsion angle is negative.
   

   angle_deg = acos(dot(b_a_orth, c_d_orth)) * 180 / pi
     % 65.609
   sign = dot(cross(b_a_orth, c_d_orth), b_c)
     % -0.91075
   if (sign < 0)
     ang_deg = -angle_deg  % -65.609
   endif
   

A related concept is the so-called dihedral angle, or more generally the angle between two planes. As long as the normal vectors to the two corresponding planes are defined, the angle between them is easy to work out.

Moreover, the method to calculate twist angles of helical nucleic acid structures in SCHNAaP and 3DNA is essentially the same.