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 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)) &gt;= 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 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".