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   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')
    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?

  • $\begingroup$ 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. $\endgroup$
    – TAR86
    Commented Jan 28, 2018 at 8:59
  • $\begingroup$ chemistry.stackexchange.com/q/55702/16683 is this relevant? $\endgroup$ Commented Jan 28, 2018 at 11:19
  • $\begingroup$ @TAR86 Thank you for your answer. However I'm trying to understand mathematical concept behind it and implement it to this script. $\endgroup$
    – Onur Ozcan
    Commented Jan 28, 2018 at 22:15
  • $\begingroup$ @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. $\endgroup$
    – Onur Ozcan
    Commented Jan 28, 2018 at 22:29

1 Answer 1


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

  • $\begingroup$ 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 $\endgroup$
    – Onur Ozcan
    Commented Jan 29, 2018 at 7:23
  • $\begingroup$ @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. $\endgroup$ Commented Jan 29, 2018 at 7:30
  • $\begingroup$ Ah I missed it. I think this might be enough, thank you again. $\endgroup$
    – Onur Ozcan
    Commented Jan 29, 2018 at 7:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.