Conversion of Z-matrix to Cartesian coordinates

To begin with, 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

```

Then, the script 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 need to perform the reverse operation and use this Z-matrix as input and define $x,y,z$ coordinates for each atom. This is converting a Z-matrix to Cartesian coordinates.


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



```
1 O          0   0   0

```

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



```
2 O          0   0   1.45335189476

```

How should I treat the 3rd and 4th atoms?


The 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 the z-coordinate of atom 2. I am not sure if I should calculate this as `z2 + distance.cos(angle)` or `z2 - distance.cos(angle)` and what it depends on, if both are possible.


For the 4th atom, I use spherical coordinates



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

```

to calculate $x,y,z$ values from the formulas



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

```

If there aren't any mistakes up to this point, how I will calculate Cartesian coordinates of the 4th atom using these $x,y,z$ values?




---


While searching I found [TMPChem's work](https://github.com/tmpchem/computational_chemistry/blob/b44ee6d2e1d2ec203811140b732cdb38f5164d76/scripts/geometry_analysis/zmat2xyz.py) on GitHub and it does exactly what I want. However, in his work, there is a mathematical part that I don't 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"?


Here is some context of how this function is used:



```
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]

```

The 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

```

Again, the only part I don't understand is `# get local axis system from 3 coordinates`. What is that function doing?




I will try to answer the question in bold from a mathematical perspective.


Basically, the first two lines in `get_local_axes` are to build the vectors $\vec r\_{12}$ and $\vec r\_{23}$. The next line builds the cross product $$\frac{\vec r\_{23}\times\vec r\_{12}}{|\vec r\_{23}\times\vec r\_{12}|},$$


The next vector is
$$\frac{\vec r\_{12}\times(\vec r\_{23}\times\vec r\_{12})}{|\vec r\_{12}\times(\vec r\_{23}\times\vec r\_{12})|},$$


This is a vector perpendicular to $\vec r\_{12}$, and lies in the plane where $\vec r\_{12}$ and $\vec r\_{23}$ lie.


The next lines are defining the three axes of the system. $\vec e\_z$ is just $\vec r\_{12}$ normalized, $\vec e\_y$ is chosen so that $yz$ is the plane where the three points (atoms) locate. and $\vec e\_x$ is chosen so that $(\vec e\_x, \vec e\_y, \vec e\_z)$ forms a orthonormal set and are the [versors](https://en.wikipedia.org/wiki/Versor_(physics)) of a right-handed coordinate system. This system is the local axis system.


So this is the meaning of "getting local axis system from 3 coordinates".


 5