What is getting local axis system from 3 coordinates?

Question

To begin with, this was originally a chemistry question before. However it changed a bit.

I wrote a script that gets Cartesian coordinates of molecule as input in the below. These are x,y,z coordinates of H2O2 molecule.

  1 O          -1.7529   -0.5188    1.3324 
  2 O          -0.4737   -0.1091    0.7774 
  3 H          -2.2902    0.1678    0.8933 
  4 H           0.0636   -0.7957    1.2164 

Than constructs a Z-matrix with them. Like this:

Z-mat :
O
O   1   1.45335189476
H   1   0.976176039452  2   96.5694760083
H   2   0.976131061897  1   96.5720363573   3   -179.995395182

Now I want to get this Z-matrix as input and define x y z coordinates for each atom. Simply converting z-matrix to Cartesian coordinates.

My question is after setting first atom as 0,0,0

1 O          0   0   0

and second one as 0,0,(distance from first) to put it on z axis

2 O          0   0   1.45335189476

How should i treat 3rd and 4th atom?

For 3rd atom must have coordinates that are something like this if read correctly:

3 H          0   distance*sin(angle)  z2+distance*cos(angle)

Taking z2 as z coordinate of atom 2.

But I am not sure if I should calculate as z2 + distance.cos(angle) or z2 - distance.cos(angle) and what it depends on if both possible.

For 4th atom should I get (r,theta,phi) like this?

r,theta,phi = (0.976,96.572,-179.995)  

to calculate x,y,z values from this spherical coordinates. Like this,

x = r * sin(theta) * cos(phi)
y = r * sin(theta) * sin(phi)
z = r * cos(theta)

If there isn't any mistake until this point;

How I will calculate Cartesian coordinates of 4th atom using these x,y,z values ?

While searching I found TMPChem's work on github. Does exactly what I want. However in his work, there is a mathematical part that I cant understand;

   # get local axis system from 3 coordinates
 def get_local_axes(coords1, coords2, coords3):
    u21 = get_u12(coords1, coords2) #calculating vector between that points 1-2
    u23 = get_u12(coords2, coords3) #calculating vector between that points 2-3
    if (abs(get_udp(u21, u23)) >= 1.0):
      print('\nError: Co-linear atoms in an internal coordinate definition')
      sys.exit()
    u23c21 = get_ucp(u23, u21) # unit cross product
    u21c23c21 = get_ucp(u21, u23c21) # unit cross product
    z = u21
    y = u21c23c21
    x = get_ucp(y, z)
    local_axes = [x, y, z]
    return local_axes

What is getting local axis system from 3 coordinates ?

Later output of this function is using here

        bond_vector = get_bond_vector(atom.rval, atom.aval, atom.tval)
        disp_vector = np.array(np.dot(bond_vector, self.atoms[i].local_axes))
        for p in range(3):
            atom.coords[p] = self.atoms[atom.rnum].coords[p] + disp_vector[p]

Bond vector definition is here:

def get_bond_vector(r, a, t):
    x = r * math.sin(a) * math.sin(t)
    y = r * math.sin(a) * math.cos(t)
    z = r * math.cos(a)
    bond_vector = [x, y, z]
    return bond_vector

The only part I dont understand is this "#get local axis system from 3 coordinates". What is the explanation of the work that function is doing there. I'm not very familiar with matrix oprerations so I'll be very grateful if the explanation is easy to understand instead of a complex and extended one.

So in short again: What is getting local axis system from 3 coordinates ?

Answer 1

First, let me say that you could have done a better job of asking this question. You told us a lot about chemistry that has nothing to do with the question of a "local axis system" and while you showed us the get_local_axis function, you didn't show us the functions it called. I had to go to github to find out what get_u12, and all those other little functions do, but you showed me get_bond_vector, which has nothing to do with the question you're asking. If you keep your questions as short and focussed as possible, you'll get a lot more responses.

Anyway, the function doesn't do anything complicated. It takes three points sat A, B, and C as input and it finds a local coordinate system of three mutually orthogonal axes, with origin at one of the points. One of the axis points in the direction of one of the other given points.

He seems to go all around the barn, though. I'm not sure why. First he constructs two unit vectors. u12 points from A to B. u23 points from B to C. u23c21 is the normalized cross product of these two vectors, so it's a unit vector perpendicular to both of them. Then u21c23c21 is the unit cross product of u21 and u23c21. Now we have three perpendicular axes. Then he does it again; I don't know why. He sets z = u21, y = u21c23c21, and computes x as the cross product of z and y. Since x is perpendicular to u21 and u21c23c21, it can only be equal to u23c21 or its negative. I don't understand the point at all.

It's just an orthogonal coordinate system with origin at one of the atoms, rotated so that one of the axes passes through one of the other two atoms.