# 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)) >= 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?

• Not really an answer, but MOLDEN converts xyz to z-matrix, IIRC, and the source should be available to you. You may have to provide a certain command line parameter to avoid a reordering of the molecule. Jan 28 '18 at 8:59
• chemistry.stackexchange.com/q/55702/16683 is this relevant? Jan 28 '18 at 11:19
• @TAR86 Thank you for your answer. However I'm trying to understand mathematical concept behind it and implement it to this script. Jan 28 '18 at 22:15
• @orthocresol I checked that question. OP asked somewhat related question trying to do a similar thing. But his question is not clear, at least I dont understand what exactly is his problem. This might be helpfull maybe. I'll try to use that source. Jan 28 '18 at 22:29

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".

• Thank you for answer. If I understand correct, dot product of this vector with bond vector(the one comes from angles) should give a scalar as displacement vector. How it has 3 component ? that allows this "for p in range(3): atom.coords[p] = self.atoms[atom.rnum].coords[p] + disp_vector[p]" to work Jan 29 '18 at 7:23
• @OnurOzcan Local axes has three components ($\vec e_x,\vec e_y,\vec e_z$), each of which is a vector, so it is in fact 3-by-3. Dot product of it with the bond vector gives the bond vector's representation (the coordinates) in the local axis system, which is composed of three elements. Jan 29 '18 at 7:30
• Ah I missed it. I think this might be enough, thank you again. Jan 29 '18 at 7:39